/*
 * case.cpp
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

#include "case_2D2t_SGA_tPrim0.h"
// MY HEADERS
#include"mros-v42.h"

namespace case_2D2t_SGA_tPrim0
{
global glob;

// *******************************************************************************************  //
// *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY ***
// *******************************************************************************************  //

#define Sqrt sqrt
#define Power pow

/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

struct parameters {
	double n;
	double J;
	double t;
	double U;
	double beta;

	unsigned int L;
};

struct variables {
	double barchiAB;
	double simchiAB;
	double DeltaS;
	double DeltaT;
	double simDeltaB;
	double mu;
	double mAF;
	double d;
    double Energy;
 };

/*
 *********************************************************
 ***************** Small useful functions ****************
 *********************************************************
*/

inline void clear(parameters &p) {
	p.n = p.J = p.t = p.U = p.beta = 0.; p.L = 0;
}

inline void clear(variables &v) {
	v.barchiAB = v.simchiAB = v.DeltaS = v.DeltaT = v.simDeltaB = v.mu = v.mAF = v.d = 0;
	v.Energy = 0;
}

inline void addTo(const variables &q, variables &v, double d=1.) {
	v.barchiAB += q.barchiAB*d;
	v.simchiAB += q.simchiAB*d;
	v.DeltaS += q.DeltaS*d;
	v.DeltaT += q.DeltaT*d;
	v.simDeltaB += q.simDeltaB*d;
	v.mu += q.mu*d;
	v.mAF += q.mAF*d;
	v.d += q.d*d;
	v.Energy += q.Energy*d;
}

inline void copyTo(variables &v, const variables &q, double d=1.) {
	v.barchiAB = q.barchiAB*d;
	v.simchiAB = q.simchiAB*d;
	v.DeltaS = q.DeltaS*d;
	v.DeltaT = q.DeltaT*d;
	v.simDeltaB = q.simDeltaB*d;
	v.mu = q.mu*d;
	v.mAF = q.mAF*d;
	v.d = q.d*d;
	v.Energy = q.Energy*d;
}

inline void copyTo(parameters &v, const parameters &q) {
	v.n = q.n;
	v.J = q.J;
	v.t = q.t;
	v.U = q.U;
	v.beta = q.beta;
	v.L = q.L;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters  p.n  p.J  p.t  p.tPrim  p.U  p.beta  p.L
	variables v; // variables  v.chiAB  v.chiS  v.chiT  v.barDelta  v.simDelta  v.mu  v.mAF  v.d
	double n, d, m, mu; // contractions

	// Abbreviations - dependend on kx, ky
	double ek;
	double etak;

public:
	double E[4];
	double dEdlbarchiAB[4];
	double dEdlsimchiAB[4];
	double dEdlDeltaS[4];
	double dEdlDeltaT[4];
	double dEdlsimDeltaB[4];
	double dEdln[4];
	double dEdlmAF[4];

	double gt, gs, dgtdm, dgtdd, dgtdn, dgsdm, dgsdd, dgsdn;
	double lbarchiAB, lsimchiAB, lDeltaS, lDeltaT, lsimDeltaB, ln, lmAF;

	Func(const parameters &p, const variables &v) {
		set(p, v);
		//refresh();
	}

	void set(const parameters &p, const variables &v) {
		copyTo(this->p, p);
		copyTo(this->v, v);
		refresh();
	}

	void set(const parameters &p) {
		copyTo(this->p, p);
		refresh();
	}

	// Gutzwiller
	double fgt();
	double fgs();
	double fgBDelta();
	double fgm();
	double fdgtdm();
	double fdgtdd();
	double fdgtdn();
	double fdgsdm();
	double fdgsdd();
	double fdgsdn();

    // Energy
    double E1(double ek, double etak);
    double E3(double ek, double etak);

    double dE1dl_barchiAB(double ek, double etak);
    double dE1dl_simchiAB(double ek, double etak);
    double dE1dl_n(double ek, double etak);
    double dE1dl_mAF(double ek, double etak);
    double dE1dl_DeltaS(double ek, double etak);
    double dE1dl_DeltaT(double ek, double etak);
    double dE1dl_simDeltaB(double ek, double etak);
    double dE1dl_mu(double ek, double etak);

    double dE3dl_barchiAB(double ek, double etak);
    double dE3dl_simchiAB(double ek, double etak);
    double dE3dl_n(double ek, double etak);
    double dE3dl_mAF(double ek, double etak);
    double dE3dl_DeltaS(double ek, double etak);
    double dE3dl_DeltaT(double ek, double etak);
    double dE3dl_simDeltaB(double ek, double etak);
    double dE3dl_mu(double ek, double etak);

	// DEBUG
	void printLambda() {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e (l)\n", lbarchiAB, lsimchiAB, lDeltaS, lDeltaT, lsimDeltaB, ln, lmAF);
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e (g)\n", gt, gs, dgtdm, dgtdd, dgtdn, dgsdm, dgsdd, dgsdn);
	}

	// REFRESH
	void refresh() {
		n = p.n;
		d = v.d;
		mu = v.mu;
		m = v.mAF;

		// obliczam funkcje Gutzwillara i poszczegolne pochodne
		gt = fgt();
		gs = fgs();
		dgtdn = fdgtdn();
		dgtdd = fdgtdd();
		dgtdm = fdgtdm();
		dgsdn = fdgsdn();
		dgsdm = fdgsdm();
		dgsdd = fdgsdd();

		// obliczanie poszczegolnych lambda:
		lbarchiAB = -gt*p.t + (3./4.)*gs*p.J*v.barchiAB;
		lsimchiAB = (1./4.)*gs*p.J*v.simchiAB;
		lDeltaS = (3./4.)*gs*p.J*v.DeltaS;
		lDeltaT = -(1./4.)*gs*p.J*v.DeltaT;
		lsimDeltaB = -(1./4.)*gs*p.J*v.simDeltaB;
		ln = -16.*p.t*v.barchiAB*dgtdn \
			- 2.*p.J*(-0.5*v.mAF*v.mAF - 3.*v.barchiAB*v.barchiAB + v.simchiAB*v.simchiAB \
                      -3.*v.DeltaS*v.DeltaS + v.DeltaT*v.DeltaT - 0.5*v.simDeltaB*v.simDeltaB)*dgsdn;
		lmAF = 2.*gs*p.J*v.mAF -16.*p.t*v.barchiAB*dgtdm \
			- 2.*p.J*(-0.5*v.mAF*v.mAF - 3.*v.barchiAB*v.barchiAB + v.simchiAB*v.simchiAB \
                      -3.*v.DeltaS*v.DeltaS + v.DeltaT*v.DeltaT - 0.5*v.simDeltaB*v.simDeltaB)*dgsdm;

		// zeby pozniej mozna bylo wypisac
		glob.lbarchiAB = lbarchiAB;
		glob.lsimchiAB = lsimchiAB;
		glob.lDeltaS = lDeltaS;
		glob.lDeltaT = lDeltaT;
		glob.lsimDeltaB = lsimDeltaB;
		glob.ln = ln;
		glob.lmAF = lmAF;
	}

	//REFRESH XY
	void refresh(double kx, double ky) {
		ek = 2.*(cos(kx) + cos(ky));
		etak = 2.*(cos(kx) - cos(ky));

		refresh();

		E[0]          = E1(ek, etak);
		dEdlbarchiAB[0]  = dE1dl_barchiAB(ek, etak);
		dEdlsimchiAB[0]  = dE1dl_simchiAB(ek, etak);
		dEdlDeltaS[0] = dE1dl_DeltaS(ek, etak);
		dEdlDeltaT[0] = dE1dl_DeltaT(ek, etak);
		dEdlsimDeltaB[0] = dE1dl_simDeltaB(ek, etak);
		dEdln[0]      = dE1dl_n(ek, etak);
		dEdlmAF[0]    = dE1dl_mAF(ek, etak);

		E[1]          = -E[0];
		dEdlbarchiAB[1]  = -dEdlbarchiAB[0];
		dEdlsimchiAB[1]  = -dEdlsimchiAB[0];
		dEdlDeltaS[1] = -dEdlDeltaS[0];
		dEdlDeltaT[1] = -dEdlDeltaT[0];
		dEdlsimDeltaB[1] = -dEdlsimDeltaB[0];
		dEdln[1]      = -dEdln[0];
		dEdlmAF[1]    = -dEdlmAF[0];

		E[2]          = E3(ek, etak);
		dEdlbarchiAB[2]  = dE3dl_barchiAB(ek, etak);
		dEdlsimchiAB[2]  = dE3dl_simchiAB(ek, etak);
		dEdlDeltaS[2] = dE3dl_DeltaS(ek, etak);
		dEdlDeltaT[2] = dE3dl_DeltaT(ek, etak);
		dEdlsimDeltaB[2] = dE3dl_simDeltaB(ek, etak);
		dEdln[2]      = dE3dl_n(ek, etak);
		dEdlmAF[2]    = dE3dl_mAF(ek, etak);

		E[3]          = -E[2];
		dEdlbarchiAB[3]  = -dEdlbarchiAB[2];
		dEdlsimchiAB[3]  = -dEdlsimchiAB[2];
		dEdlDeltaS[3] = -dEdlDeltaS[2];
		dEdlDeltaT[3] = -dEdlDeltaT[2];
		dEdlsimDeltaB[3] = -dEdlsimDeltaB[2];
		dEdln[3]      = -dEdln[2];
		dEdlmAF[3]    = -dEdlmAF[2];
	}
};


/////////////////
///////////////// FullSimplify
/////////////////

double Func::fgt() {
	return (2*(2*Power(d,2) - n)*(Sqrt(((2 + m - n)*(-1 - Power(d,2) + n))/(-2 + m + n)) +
		       Sqrt((Power(d,2)*(-m + n))/(m + n)))*
		     (Sqrt(((-1 - Power(d,2) + n)*(-2 + m + n))/(2 + m - n)) + Sqrt((Power(d,2)*(m + n))/(-m + n))))/
		   (-Power(m,2) + (-2 + n)*n);
}

double Func::fgs() {
	return (4*Power(-2*Power(d,2) + n,2))/Power(Power(m,2) - (-2 + n)*n,2);
}

double Func::fgBDelta() {
    //n = 1.-p.sdelta;
    double r = n/2. + v.mAF/2.;
    double w = n/2. - v.mAF/2.;
    double d=v.d*v.d;
	double part1 = (n-2.*d)/(2.*(n-2.*r*w));
	double part2 = sqrt((1.-w)*(1.-n+d)/(1.-r)) + sqrt(w*d/r);
	double part3 = sqrt((1.-r)*(1.-n+d)/(1.-w)) + sqrt(r*d/w);
	return part1*(part2*part2 + part3*part3);
}

double Func::fgm() {
   double r = n/2. + v.mAF/2.;
   double w = n/2. - v.mAF/2.;
   double d = v.d*v.d;
   return (n-2.*d)/(n-2.*w*r);
}

double Func::fdgsdm() {
	return (-16*m*Power(-2*Power(d,2) + n,2))/Power(Power(m,2) - (-2 + n)*n,3);
}

double Func::fdgsdd() {
	return (64*Power(d,3) - 32*d*n)/Power(Power(m,2) - (-2 + n)*n,2);
}

double Func::fdgsdn() {
	return (8*(-2*Power(d,2) + n)*(Power(m,2) - 4*Power(d,2)*(-1 + n) + Power(n,2)))/Power(Power(m,2) - (-2 + n)*n,3);
}

double Func::fdgtdm() {
	return (4*(-2*Power(d,2) + n)*(-(m*(Sqrt(((2 + m - n)*(-1 - Power(d,2) + n))/(-2 + m + n)) +
            Sqrt((Power(d,2)*(-m + n))/(m + n)))*
          (Sqrt(((-1 - Power(d,2) + n)*(-2 + m + n))/(2 + m - n)) + Sqrt((Power(d,2)*(m + n))/(-m + n)))) +
       (Power(m,4)*(Sqrt(((-1 - Power(d,2) + n)*(-2 + m + n))/(2 + m - n))*
             Sqrt((Power(d,2)*(-m + n))/(m + n)) -
            Sqrt(((2 + m - n)*(-1 - Power(d,2) + n))/(-2 + m + n))*Sqrt((Power(d,2)*(m + n))/(-m + n))))/
        ((m - n)*(2 + m - n)*(-2 + m + n)*(m + n)) -
       (Power(-2 + n,2)*Power(n,2)*(Sqrt(((-1 - Power(d,2) + n)*(-2 + m + n))/(2 + m - n))*
             Sqrt((Power(d,2)*(-m + n))/(m + n)) -
            Sqrt(((2 + m - n)*(-1 - Power(d,2) + n))/(-2 + m + n))*Sqrt((Power(d,2)*(m + n))/(-m + n))))/
        ((m - n)*(2 + m - n)*(-2 + m + n)*(m + n))))/Power(Power(m,2) - (-2 + n)*n,2);
}

double Func::fdgtdd() {
	return (2*(-4*d*(Sqrt(((2 + m - n)*(-1 - Power(d,2) + n))/(-2 + m + n)) + Sqrt((Power(d,2)*(-m + n))/(m + n)))*
	        (Sqrt(((-1 - Power(d,2) + n)*(-2 + m + n))/(2 + m - n)) + Sqrt((Power(d,2)*(m + n))/(-m + n))) +
	       (-2*Power(d,2) + n)*((d*Sqrt(((2 + m - n)*(-1 - Power(d,2) + n))/(-2 + m + n)))/
	           (1 + Power(d,2) - n) + Sqrt((Power(d,2)*(-m + n))/(m + n))/d)*
	        (Sqrt(((-1 - Power(d,2) + n)*(-2 + m + n))/(2 + m - n)) + Sqrt((Power(d,2)*(m + n))/(-m + n))) +
	       (-2*Power(d,2) + n)*(Sqrt(((2 + m - n)*(-1 - Power(d,2) + n))/(-2 + m + n)) +
	          Sqrt((Power(d,2)*(-m + n))/(m + n)))*
	        ((d*Sqrt(((-1 - Power(d,2) + n)*(-2 + m + n))/(2 + m - n)))/(1 + Power(d,2) - n) +
	          Sqrt((Power(d,2)*(m + n))/(-m + n))/d)))/(Power(m,2) - (-2 + n)*n);
}

double Func::fdgtdn() {
	return (4*(((Power(m,2) - (-2 + n)*n)*(Sqrt(((2 + m - n)*(-1 - Power(d,2) + n))/(-2 + m + n)) +
            Sqrt((Power(d,2)*(-m + n))/(m + n)))*
          (Sqrt(((-1 - Power(d,2) + n)*(-2 + m + n))/(2 + m - n)) + Sqrt(-((Power(d,2)*(m + n))/(m - n)))))/
        2. - (1 - n)*(-2*Power(d,2) + n)*(Sqrt(((2 + m - n)*(-1 - Power(d,2) + n))/(-2 + m + n)) +
          Sqrt((Power(d,2)*(-m + n))/(m + n)))*
        (Sqrt(((-1 - Power(d,2) + n)*(-2 + m + n))/(2 + m - n)) + Sqrt((Power(d,2)*(m + n))/(-m + n))) +
       ((-2*Power(d,2) + n)*(Power(m,2) - (-2 + n)*n)*
          ((Power(m,2) + 2*m*(1 + Power(d,2) - n) - Power(-2 + n,2))/
             (Sqrt(((2 + m - n)*(-1 - Power(d,2) + n))/(-2 + m + n))*Power(-2 + m + n,2)) +
            (2*Power(d,2)*m)/(Sqrt((Power(d,2)*(-m + n))/(m + n))*Power(m + n,2)))*
          (Sqrt(((-1 - Power(d,2) + n)*(-2 + m + n))/(2 + m - n)) + Sqrt((Power(d,2)*(m + n))/(-m + n))))/4.\
        + ((-2*Power(d,2) + n)*(Power(m,2) - (-2 + n)*n)*
          (Sqrt(((2 + m - n)*(-1 - Power(d,2) + n))/(-2 + m + n)) + Sqrt((Power(d,2)*(-m + n))/(m + n)))*
          ((Power(m,2) - 2*m*(1 + Power(d,2) - n) - Power(-2 + n,2))/
             (Power(2 + m - n,2)*Sqrt(((-1 - Power(d,2) + n)*(-2 + m + n))/(2 + m - n))) +
            (2*m*Sqrt(-((Power(d,2)*(m + n))/(m - n))))/(Power(m,2) - Power(n,2))))/4.))/
   Power(Power(m,2) - (-2 + n)*n,2);
}


double Func::E1(double ek, double etak) {
	return -Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		      2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2) -
		      2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		         4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		          (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2))))/2.;
}

double Func::E3(double ek, double etak) {
	return -Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		      2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2) +
		      2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		         4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		          (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2))))/2.;
}

double Func::dE1dl_barchiAB(double ek, double etak) {
	return (2*Power(ek,2)*lbarchiAB*(Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2) -
		       Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		         4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		          (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
		   (Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		       4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		        (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
		     Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
		       Power(ln + 2*mu,2) - 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))));
}

double Func::dE3dl_barchiAB(double ek, double etak) {
	return (-2*Power(ek,2)*lbarchiAB*(Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2) +
		       Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		         4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		          (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
		   (Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		       4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		        (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
		     Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
		       Power(ln + 2*mu,2) + 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))));
}

double Func::dE1dl_simchiAB(double ek, double etak) {
	return (2*Power(ek,2)*lsimchiAB*(Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2) -
		       Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		         4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		          (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
		   (Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		       4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		        (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
		     Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
		       Power(ln + 2*mu,2) - 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))));
}

double Func::dE3dl_simchiAB(double ek, double etak) {
	return (-2*Power(ek,2)*lsimchiAB*(Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2) +
		       Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		         4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		          (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
		   (Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		       4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		        (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
		     Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
		       Power(ln + 2*mu,2) + 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))));
}

double Func::dE1dl_n(double ek, double etak) {
	return (Power(etak,2)*lmAF*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) -
		     (ln + 2*mu)*(-Power(lmAF,2) - 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		        Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		          4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
		   (2.*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		       4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		        (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
		     Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
		       Power(ln + 2*mu,2) - 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))));
}

double Func::dE3dl_n(double ek, double etak) {
	return (Power(etak,2)*lmAF*(-4*lDeltaS*lDeltaT + Power(lsimDeltaB,2)) -
		     (ln + 2*mu)*(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		        Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		          4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
		   (2.*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		       4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		        (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
		     Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
		       Power(ln + 2*mu,2) + 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))));
}

double Func::dE1dl_mAF(double ek, double etak) {
	return (Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2))*(ln + 2*mu) +
		     lmAF*(Power(ln + 2*mu,2) - Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
		   (2.*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		       4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		        (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
		     Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
		       Power(ln + 2*mu,2) - 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))));
}

double Func::dE3dl_mAF(double ek, double etak) {
	return (-(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2))*(ln + 2*mu)) - lmAF*Power(ln + 2*mu,2) -
		     lmAF*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		        4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		         (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2))))/
		   (2.*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		       4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		        (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
		     Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
		       Power(ln + 2*mu,2) + 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))));
}

double Func::dE1dl_DeltaS(double ek, double etak) {
	return (2*Power(etak,2)*(Power(etak,2)*lDeltaT*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		       lDeltaT*lmAF*(ln + 2*mu) - lDeltaS*Sqrt(Power(Power(etak,2)*
		             (4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		          4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
		   (Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		       4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		        (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
		     Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
		       Power(ln + 2*mu,2) - 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))));
}

double Func::dE3dl_DeltaS(double ek, double etak) {
	return (-2*Power(etak,2)*(Power(etak,2)*lDeltaT*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		       lDeltaT*lmAF*(ln + 2*mu) + lDeltaS*Sqrt(Power(Power(etak,2)*
		             (4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		          4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
		   (Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		       4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		        (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
		     Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
		       Power(ln + 2*mu,2) + 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))));
}

double Func::dE1dl_DeltaT(double ek, double etak) {
	return (2*Power(etak,2)*(4*Power(ek,2)*lDeltaT*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       Power(etak,2)*lDeltaS*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lDeltaS*lmAF*(ln + 2*mu) -
		       lDeltaT*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		          4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
		   (Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		       4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		        (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
		     Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
		       Power(ln + 2*mu,2) - 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))));
}

double Func::dE3dl_DeltaT(double ek, double etak) {
	return (-2*Power(etak,2)*(4*Power(ek,2)*lDeltaT*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       Power(etak,2)*lDeltaS*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lDeltaS*lmAF*(ln + 2*mu) +
		       lDeltaT*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		          4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
		   (Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
		       4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		        (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
		     Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
		       2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
		       Power(ln + 2*mu,2) + 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
		            lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
		           (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))));
}

double Func::dE1dl_simDeltaB(double ek, double etak) {
	return -((Power(etak,2)*lsimDeltaB*(-2*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
	         Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu) +
	         Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
	           4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
	            (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
	     (Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
	         4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
	          (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
	       Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
	         2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
	         Power(ln + 2*mu,2) - 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
	              lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
	             (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2))))));
}

double Func::dE3dl_simDeltaB(double ek, double etak) {
	return -((Power(etak,2)*lsimDeltaB*(2*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
	         Power(etak,2)*(-4*lDeltaS*lDeltaT + Power(lsimDeltaB,2)) - lmAF*(ln + 2*mu) +
	         Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
	           4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
	            (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))))/
	     (Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) + lmAF*(ln + 2*mu),2) +
	         4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
	          (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2)))*
	       Sqrt(Power(lmAF,2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2)) +
	         2*Power(etak,2)*(2*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lsimDeltaB,2)) +
	         Power(ln + 2*mu,2) + 2*Sqrt(Power(Power(etak,2)*(4*lDeltaS*lDeltaT - Power(lsimDeltaB,2)) +
	              lmAF*(ln + 2*mu),2) + 4*Power(ek,2)*(Power(lbarchiAB,2) + Power(lsimchiAB,2))*
	             (Power(etak,2)*(4*Power(lDeltaT,2) + Power(lsimDeltaB,2)) + Power(ln + 2*mu,2))))));
}

/*
double Func::fgt() {
	return ((Sqrt(((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(1 - m/2. - n/2.)) + Sqrt((Power(d,2)*(-m/2. + n/2.))/(m/2. + n/2.)))*
		     (Sqrt(((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(1 + m/2. - n/2.)) + Sqrt((Power(d,2)*(m/2. + n/2.))/(-m/2. + n/2.)))*(-2*Power(d,2) + n))/(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n);
}

double Func::fgs() {
	return Power(-2*Power(d,2) + n,2)/Power(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n,2);
}

double Func::fgBDelta() {
    //n = 1.-p.sdelta;
    double r = n/2. + v.mAF/2.;
    double w = n/2. - v.mAF/2.;
	double part1 = (n-2.*v.d)/(2.*(n-2.*r*w));
	double part2 = sqrt((1.-w)*(1.-n+v.d)/(1.-r)) + sqrt(w*v.d/r);
	double part3 = sqrt((1.-r)*(1.-n+v.d)/(1.-w)) + sqrt(r*v.d/w);
	return part1*(part2*part2 + part3*part3);
}

double Func::fgm() {
   double r = n/2. + v.mAF/2.;
   double w = n/2. - v.mAF/2.;
   return (n-2.*v.d)/(n-2.*w*r);
}

double Func::fdgsdm() {
	return (-2*m*Power(-2*Power(d,2) + n,2))/Power(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n,3);
}

double Func::fdgsdd() {
	return (-8*d*(-2*Power(d,2) + n))/Power(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n,2);
}

double Func::fdgsdn() {
	return (-2*(1 - n)*Power(-2*Power(d,2) + n,2))/Power(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n,3) + (2*(-2*Power(d,2) + n))/Power(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n,2);
}

double Func::fdgtdm() {
	return -((m*(Sqrt(((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(1 - m/2. - n/2.)) + Sqrt((Power(d,2)*(-m/2. + n/2.))/(m/2. + n/2.)))*
	        (Sqrt(((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(1 + m/2. - n/2.)) + Sqrt((Power(d,2)*(m/2. + n/2.))/(-m/2. + n/2.)))*(-2*Power(d,2) + n))/Power(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n,2)) +
	   ((Sqrt(((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(1 - m/2. - n/2.)) + Sqrt((Power(d,2)*(-m/2. + n/2.))/(m/2. + n/2.)))*
	      ((-((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(2.*Power(1 + m/2. - n/2.,2)) - (1 + Power(d,2) - n)/(2.*(1 + m/2. - n/2.)))/(2.*Sqrt(((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(1 + m/2. - n/2.))) +
	        (Power(d,2)/(2.*(-m/2. + n/2.)) + (Power(d,2)*(m/2. + n/2.))/(2.*Power(-m/2. + n/2.,2)))/(2.*Sqrt((Power(d,2)*(m/2. + n/2.))/(-m/2. + n/2.))))*(-2*Power(d,2) + n))/(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n) +
	   ((((1 + Power(d,2) - n)/(2.*(1 - m/2. - n/2.)) + ((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(2.*Power(1 - m/2. - n/2.,2)))/(2.*Sqrt(((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(1 - m/2. - n/2.))) +
	        (-(Power(d,2)*(-m/2. + n/2.))/(2.*Power(m/2. + n/2.,2)) - Power(d,2)/(2.*(m/2. + n/2.)))/(2.*Sqrt((Power(d,2)*(-m/2. + n/2.))/(m/2. + n/2.))))*
	      (Sqrt(((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(1 + m/2. - n/2.)) + Sqrt((Power(d,2)*(m/2. + n/2.))/(-m/2. + n/2.)))*(-2*Power(d,2) + n))/(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n);
}

double Func::fdgtdd() {
	return (-4*d*(Sqrt(((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(1 - m/2. - n/2.)) + Sqrt((Power(d,2)*(-m/2. + n/2.))/(m/2. + n/2.)))*
		      (Sqrt(((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(1 + m/2. - n/2.)) + Sqrt((Power(d,2)*(m/2. + n/2.))/(-m/2. + n/2.))))/(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n) +
		   ((Sqrt(((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(1 - m/2. - n/2.)) + Sqrt((Power(d,2)*(-m/2. + n/2.))/(m/2. + n/2.)))*
		      ((d*(1 - m/2. - n/2.))/(Sqrt(((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(1 + m/2. - n/2.))*(1 + m/2. - n/2.)) + (d*(m/2. + n/2.))/((-m/2. + n/2.)*Sqrt((Power(d,2)*(m/2. + n/2.))/(-m/2. + n/2.))))*(-2*Power(d,2) + n))/
		    (-2*(-m/2. + n/2.)*(m/2. + n/2.) + n) + (((d*(1 + m/2. - n/2.))/((1 - m/2. - n/2.)*Sqrt(((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(1 - m/2. - n/2.))) +
		        (d*(-m/2. + n/2.))/(Sqrt((Power(d,2)*(-m/2. + n/2.))/(m/2. + n/2.))*(m/2. + n/2.)))*(Sqrt(((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(1 + m/2. - n/2.)) + Sqrt((Power(d,2)*(m/2. + n/2.))/(-m/2. + n/2.)))*
		      (-2*Power(d,2) + n))/(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n);
}

double Func::fdgtdn() {
	return -(((Sqrt(((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(1 - m/2. - n/2.)) + Sqrt((Power(d,2)*(-m/2. + n/2.))/(m/2. + n/2.)))*
	        (Sqrt(((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(1 + m/2. - n/2.)) + Sqrt((Power(d,2)*(m/2. + n/2.))/(-m/2. + n/2.)))*(1 - n)*(-2*Power(d,2) + n))/Power(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n,2)) +
	   ((Sqrt(((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(1 - m/2. - n/2.)) + Sqrt((Power(d,2)*(-m/2. + n/2.))/(m/2. + n/2.)))*
	      (Sqrt(((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(1 + m/2. - n/2.)) + Sqrt((Power(d,2)*(m/2. + n/2.))/(-m/2. + n/2.))))/(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n) +
	   ((Sqrt(((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(1 - m/2. - n/2.)) + Sqrt((Power(d,2)*(-m/2. + n/2.))/(m/2. + n/2.)))*
	      ((((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(2.*Power(1 + m/2. - n/2.,2)) - (1 + Power(d,2) - n)/(2.*(1 + m/2. - n/2.)) - (1 - m/2. - n/2.)/(1 + m/2. - n/2.))/
	         (2.*Sqrt(((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(1 + m/2. - n/2.))) + (Power(d,2)/(2.*(-m/2. + n/2.)) - (Power(d,2)*(m/2. + n/2.))/(2.*Power(-m/2. + n/2.,2)))/(2.*Sqrt((Power(d,2)*(m/2. + n/2.))/(-m/2. + n/2.))))*
	      (-2*Power(d,2) + n))/(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n) + (((-(1 + Power(d,2) - n)/(2.*(1 - m/2. - n/2.)) + ((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(2.*Power(1 - m/2. - n/2.,2)) -
	           (1 + m/2. - n/2.)/(1 - m/2. - n/2.))/(2.*Sqrt(((1 + Power(d,2) - n)*(1 + m/2. - n/2.))/(1 - m/2. - n/2.))) +
	        (-(Power(d,2)*(-m/2. + n/2.))/(2.*Power(m/2. + n/2.,2)) + Power(d,2)/(2.*(m/2. + n/2.)))/(2.*Sqrt((Power(d,2)*(-m/2. + n/2.))/(m/2. + n/2.))))*
	      (Sqrt(((1 + Power(d,2) - n)*(1 - m/2. - n/2.))/(1 + m/2. - n/2.)) + Sqrt((Power(d,2)*(m/2. + n/2.))/(-m/2. + n/2.)))*(-2*Power(d,2) + n))/(-2*(-m/2. + n/2.)*(m/2. + n/2.) + n);
}


double Func::E1(double ek, double etak) {
	return -Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
		      4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
		      2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) -
		      2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
		         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
		         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
		         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
		         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
		         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
		         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
		         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
		         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
		         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
		         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
		         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
		         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/2.;
}

double Func::E3(double ek, double etak) {
	return -Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
		      4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
		      2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) +
		      2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
		         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
		         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
		         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
		         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
		         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
		         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
		         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
		         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
		         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
		         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
		         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
		         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/2.;
}

double Func::dE1dl_barchiAB(double ek, double etak) {
	return -(8*Power(ek,2)*lbarchiAB - (32*Power(ek,2)*Power(etak,2)*lbarchiAB*Power(lDeltaT,2) +
	         8*Power(ek,2)*lbarchiAB*Power(ln,2) + 8*Power(ek,2)*Power(etak,2)*lbarchiAB*Power(lsimDeltaB,2) +
	         32*Power(ek,2)*lbarchiAB*ln*mu + 32*Power(ek,2)*lbarchiAB*Power(mu,2))/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) -
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE3dl_barchiAB(double ek, double etak) {
	return -(8*Power(ek,2)*lbarchiAB + (32*Power(ek,2)*Power(etak,2)*lbarchiAB*Power(lDeltaT,2) +
	         8*Power(ek,2)*lbarchiAB*Power(ln,2) + 8*Power(ek,2)*Power(etak,2)*lbarchiAB*Power(lsimDeltaB,2) +
	         32*Power(ek,2)*lbarchiAB*ln*mu + 32*Power(ek,2)*lbarchiAB*Power(mu,2))/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) +
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE1dl_simchiAB(double ek, double etak) {
	return -(8*Power(ek,2)*lsimchiAB - (32*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*lsimchiAB +
	         8*Power(ek,2)*Power(ln,2)*lsimchiAB + 8*Power(ek,2)*Power(etak,2)*lsimchiAB*Power(lsimDeltaB,2) +
	         32*Power(ek,2)*ln*lsimchiAB*mu + 32*Power(ek,2)*lsimchiAB*Power(mu,2))/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) -
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE3dl_simchiAB(double ek, double etak) {
	return -(8*Power(ek,2)*lsimchiAB + (32*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*lsimchiAB +
	         8*Power(ek,2)*Power(ln,2)*lsimchiAB + 8*Power(ek,2)*Power(etak,2)*lsimchiAB*Power(lsimDeltaB,2) +
	         32*Power(ek,2)*ln*lsimchiAB*mu + 32*Power(ek,2)*lsimchiAB*Power(mu,2))/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) +
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE1dl_n(double ek, double etak) {
	return -(2*ln + 4*mu - (8*Power(etak,2)*lDeltaS*lDeltaT*lmAF + 8*Power(ek,2)*Power(lbarchiAB,2)*ln +
	         2*Power(lmAF,2)*ln + 8*Power(ek,2)*ln*Power(lsimchiAB,2) -
	         2*Power(etak,2)*lmAF*Power(lsimDeltaB,2) + 16*Power(ek,2)*Power(lbarchiAB,2)*mu +
	         4*Power(lmAF,2)*mu + 16*Power(ek,2)*Power(lsimchiAB,2)*mu)/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) -
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE3dl_n(double ek, double etak) {
	return -(2*ln + 4*mu + (8*Power(etak,2)*lDeltaS*lDeltaT*lmAF + 8*Power(ek,2)*Power(lbarchiAB,2)*ln +
	         2*Power(lmAF,2)*ln + 8*Power(ek,2)*ln*Power(lsimchiAB,2) -
	         2*Power(etak,2)*lmAF*Power(lsimDeltaB,2) + 16*Power(ek,2)*Power(lbarchiAB,2)*mu +
	         4*Power(lmAF,2)*mu + 16*Power(ek,2)*Power(lsimchiAB,2)*mu)/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) +
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE1dl_mAF(double ek, double etak) {
	return -(2*lmAF - (8*Power(etak,2)*lDeltaS*lDeltaT*ln + 2*lmAF*Power(ln,2) -
	         2*Power(etak,2)*ln*Power(lsimDeltaB,2) + 16*Power(etak,2)*lDeltaS*lDeltaT*mu + 8*lmAF*ln*mu -
	         4*Power(etak,2)*Power(lsimDeltaB,2)*mu + 8*lmAF*Power(mu,2))/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) -
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE3dl_mAF(double ek, double etak) {
	return -(2*lmAF + (8*Power(etak,2)*lDeltaS*lDeltaT*ln + 2*lmAF*Power(ln,2) -
	         2*Power(etak,2)*ln*Power(lsimDeltaB,2) + 16*Power(etak,2)*lDeltaS*lDeltaT*mu + 8*lmAF*ln*mu -
	         4*Power(etak,2)*Power(lsimDeltaB,2)*mu + 8*lmAF*Power(mu,2))/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) +
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE1dl_DeltaS(double ek, double etak) {
	return -(8*Power(etak,2)*lDeltaS - (32*Power(etak,4)*lDeltaS*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaT*lmAF*ln -
	         8*Power(etak,4)*lDeltaT*Power(lsimDeltaB,2) + 16*Power(etak,2)*lDeltaT*lmAF*mu)/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) -
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE3dl_DeltaS(double ek, double etak) {
	return -(8*Power(etak,2)*lDeltaS + (32*Power(etak,4)*lDeltaS*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaT*lmAF*ln -
	         8*Power(etak,4)*lDeltaT*Power(lsimDeltaB,2) + 16*Power(etak,2)*lDeltaT*lmAF*mu)/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) +
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE1dl_DeltaT(double ek, double etak) {
	return -(8*Power(etak,2)*lDeltaT - (32*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*lDeltaT +
	         32*Power(etak,4)*Power(lDeltaS,2)*lDeltaT + 8*Power(etak,2)*lDeltaS*lmAF*ln +
	         32*Power(ek,2)*Power(etak,2)*lDeltaT*Power(lsimchiAB,2) -
	         8*Power(etak,4)*lDeltaS*Power(lsimDeltaB,2) + 16*Power(etak,2)*lDeltaS*lmAF*mu)/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) -
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE3dl_DeltaT(double ek, double etak) {
	return -(8*Power(etak,2)*lDeltaT + (32*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*lDeltaT +
	         32*Power(etak,4)*Power(lDeltaS,2)*lDeltaT + 8*Power(etak,2)*lDeltaS*lmAF*ln +
	         32*Power(ek,2)*Power(etak,2)*lDeltaT*Power(lsimchiAB,2) -
	         8*Power(etak,4)*lDeltaS*Power(lsimDeltaB,2) + 16*Power(etak,2)*lDeltaS*lmAF*mu)/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) +
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE1dl_simDeltaB(double ek, double etak) {
	return -(4*Power(etak,2)*lsimDeltaB - (8*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*lsimDeltaB -
	         16*Power(etak,4)*lDeltaS*lDeltaT*lsimDeltaB - 4*Power(etak,2)*lmAF*ln*lsimDeltaB +
	         8*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*lsimDeltaB + 4*Power(etak,4)*Power(lsimDeltaB,3) -
	         8*Power(etak,2)*lmAF*lsimDeltaB*mu)/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) -
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}

double Func::dE3dl_simDeltaB(double ek, double etak) {
	return -(4*Power(etak,2)*lsimDeltaB + (8*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*lsimDeltaB -
	         16*Power(etak,4)*lDeltaS*lDeltaT*lsimDeltaB - 4*Power(etak,2)*lmAF*ln*lsimDeltaB +
	         8*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*lsimDeltaB + 4*Power(etak,4)*Power(lsimDeltaB,3) -
	         8*Power(etak,2)*lmAF*lsimDeltaB*mu)/
	       Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	         16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	         4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	         16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	         8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) - 2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	         4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	         Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	         16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	         16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	         16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2)))/
	   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lDeltaS,2) +
	       4*Power(etak,2)*Power(lDeltaT,2) + Power(lmAF,2) + Power(ln,2) + 4*Power(ek,2)*Power(lsimchiAB,2) +
	       2*Power(etak,2)*Power(lsimDeltaB,2) + 4*ln*mu + 4*Power(mu,2) +
	       2*Sqrt(16*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lDeltaT,2) +
	          16*Power(etak,4)*Power(lDeltaS,2)*Power(lDeltaT,2) + 8*Power(etak,2)*lDeltaS*lDeltaT*lmAF*ln +
	          4*Power(ek,2)*Power(lbarchiAB,2)*Power(ln,2) + Power(lmAF,2)*Power(ln,2) +
	          16*Power(ek,2)*Power(etak,2)*Power(lDeltaT,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(ln,2)*Power(lsimchiAB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lbarchiAB,2)*Power(lsimDeltaB,2) -
	          8*Power(etak,4)*lDeltaS*lDeltaT*Power(lsimDeltaB,2) -
	          2*Power(etak,2)*lmAF*ln*Power(lsimDeltaB,2) +
	          4*Power(ek,2)*Power(etak,2)*Power(lsimchiAB,2)*Power(lsimDeltaB,2) +
	          Power(etak,4)*Power(lsimDeltaB,4) + 16*Power(etak,2)*lDeltaS*lDeltaT*lmAF*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*ln*mu + 4*Power(lmAF,2)*ln*mu +
	          16*Power(ek,2)*ln*Power(lsimchiAB,2)*mu - 4*Power(etak,2)*lmAF*Power(lsimDeltaB,2)*mu +
	          16*Power(ek,2)*Power(lbarchiAB,2)*Power(mu,2) + 4*Power(lmAF,2)*Power(mu,2) +
	          16*Power(ek,2)*Power(lsimchiAB,2)*Power(mu,2))));
}
*/

/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	void onePoint(variables &y, int i, int j);

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->set(p);
	}

	void run(variables &w);
};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	f->refresh(x, y);
	if(DEEP_DEBUG) { f->printLambda(); }

	double Etan[4] = { (1-tanh(p.beta*f->E[0]/2.))/2., (1-tanh(p.beta*f->E[1]/2.))/2., (1-tanh(p.beta*f->E[2]/2.))/2., (1-tanh(p.beta*f->E[3]/2.))/2. };

	for(int i=0; i<4; i++) {
		w.barchiAB += Etan[i]*f->dEdlbarchiAB[i];
		w.simchiAB += Etan[i]*f->dEdlsimchiAB[i];
		w.DeltaS += Etan[i]*f->dEdlDeltaS[i];
		w.DeltaT += Etan[i]*f->dEdlDeltaT[i];
		w.simDeltaB += Etan[i]*f->dEdlsimDeltaB[i];
		w.mu += Etan[i]*f->dEdln[i];
		w.mAF += Etan[i]*f->dEdlmAF[i];
		// w.d w tym miejscu jeszcze nie

        //if(f->E[i]>0) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
        //if(f->E[i]<0) { w.Energy += f->E[i]; }

        // it would be safer have following code - warning: 30 is chosen arbitrarily
		if(p.beta*f->E[i] > 30.) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
		else if(p.beta*f->E[i] < -30) { w.Energy += f->E[i]; }
		else { w.Energy += -log(1 + exp(-p.beta*f->E[i]))/p.beta; }
	}
}

void Equa::run(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation
	unsigned int max = p.L/2; // attention: L should be even!

	// SERIAL CASE
	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
	// * edges (2 times)
	for(unsigned int i=1; i<max; i++) {
		onePoint(temp, i, 0);
		onePoint(temp, 0, i);
		onePoint(temp, i, max-i);
	}
   addTo(temp, w, 2);
   clear(temp);
   // * interior area (4 times) - SERIAL CASE
   for(unsigned int i=1; i<max; i++) {
		for(unsigned int j=1; j<i; j++) {
			onePoint(temp, max-i, j);
     	}
	}
	addTo(temp, w, 4);
	clear(temp);
	// * missing corners
	onePoint(w, max, 0);
	//onePoint(w, 0, max);
	onePoint(w, 0, 0);

	// what was out of sumation
	double L2 = (double)(p.L*p.L);
	w.barchiAB /= L2;
	w.simchiAB /= L2;
	w.DeltaS /= L2;
	w.DeltaT /= L2;
	w.simDeltaB /= L2;
	w.mu /= L2;
	w.mAF /= L2;
    w.Energy /= L2;

	w.barchiAB += 8.*v.barchiAB;
	w.simchiAB += -8.*v.simchiAB;
	w.DeltaS += 8.*v.DeltaS;
	w.DeltaT += 8.*v.DeltaT;
	w.simDeltaB += -4.*v.simDeltaB;
	w.mu += -0.5*(1-p.n);
	w.mAF += 0.5*v.mAF;

	// calculate dd
	w.d = 2.*v.d*p.U + 8.*p.t*v.barchiAB*f->dgtdd \
			+ p.J*(-0.5*v.mAF*v.mAF - 3.*v.barchiAB*v.barchiAB + v.simchiAB*v.simchiAB \
                    -3.*v.DeltaS*v.DeltaS + v.DeltaT*v.DeltaT - 0.5*v.simDeltaB*v.simDeltaB)*f->dgsdd;

    v.Energy = w.Energy \
    		+ 8.*f->gt*p.t*v.barchiAB \
    		+ f->gs*p.J*(-0.5*v.mAF*v.mAF - 3.*v.barchiAB*v.barchiAB + v.simchiAB*v.simchiAB \
                    -3.*v.DeltaS*v.DeltaS + v.DeltaT*v.DeltaT - 0.5*v.simDeltaB*v.simDeltaB) \
    		+ 8.*(f->lbarchiAB*v.barchiAB - f->lsimchiAB*v.simchiAB + f->lDeltaS*v.DeltaS + f->lDeltaT*v.DeltaT - 0.5*f->lsimDeltaB*v.simDeltaB) \
    		+ p.U*v.d*v.d + 0.5*f->ln*(p.n-1) + 0.5*f->lmAF*v.mAF -v.mu;

	//if(w.chiAB != w.chiAB) { w.chiAB = 1000; }
	//if(w.chiS != w.chiS) { w.chiS = 1000; }
	//if(w.chiT != w.chiT) { w.chiT = 1000; }
	//if(w.barDelta != w.barDelta) { w.barDelta = 1000; }
	//if(w.simDelta != w.simDelta) { w.simDelta = 1000; }
	//if(w.mu != w.mu) { w.mu = 1000; }
	//if(w.mAF != w.mAF) { w.mAF = 1000; }

    // results
    case_2D2t_SGA_tPrim0::glob.Energy = v.Energy + v.mu*p.n;
	glob.gt = f->fgt();
	glob.gs = f->fgs();
	glob.gBDelta = f->fgBDelta();
	glob.gm = f->fgm();
}

// ***********************************************************************************  //
// *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY ***
// ***********************************************************************************  //


int function(const gsl_vector * x, void *params, gsl_vector * f) {

	struct parameters p;
	p.n = ((Par *) params)->n();
	p.J = ((Par *) params)->J();
	p.t = ((Par *) params)->t();
	p.U = ((Par *) params)->U();
	p.beta = ((Par *) params)->beta();
	p.L = ((Par *) params)->L();

	variables v;
	v.barchiAB = gsl_vector_get (x, 0);
	v.simchiAB = gsl_vector_get (x, 1);
	v.DeltaS = gsl_vector_get (x, 2);
	v.DeltaT = gsl_vector_get (x, 3);
	v.simDeltaB = gsl_vector_get (x, 4);
	v.mu = gsl_vector_get (x, 5);
	v.mAF = gsl_vector_get (x, 6);
	v.d = gsl_vector_get (x, 7);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15ee\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5),
				gsl_vector_get (x, 6), gsl_vector_get (x, 7));
	}

    // It's start of the calculations

    Equa E(p, v);
    variables w;
    E.run(w);

    // It's end of the calculations

	gsl_vector_set (f, 0, w.barchiAB);
	gsl_vector_set (f, 1, w.simchiAB);
	gsl_vector_set (f, 2, w.DeltaS);
	gsl_vector_set (f, 3, w.DeltaT);
	gsl_vector_set (f, 4, w.simDeltaB);
	gsl_vector_set (f, 5, w.mu);
	gsl_vector_set (f, 6, w.mAF);
	gsl_vector_set (f, 7, w.d);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5),
				gsl_vector_get (x, 6), gsl_vector_get (x, 7));
		printf("###  %24.15e# %25.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e | E = %24.15e\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2),
			gsl_vector_get (f, 3), gsl_vector_get (f, 4), gsl_vector_get (f, 5),
			gsl_vector_get (f, 6), gsl_vector_get (f, 7), case_2D2t_SGA_tPrim0::glob.Energy);
		printf("\n");
	}

    return GSL_SUCCESS;
}

// ************* CLASSES ************* //

// Par

void Par::print(FILE *file) const {
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
}

void Par::print_verboten(FILE *file) const {
	fprintf(file, "PARAMETERS: n=%f, J=%f, t=%f, U=%f, beta=%f, L=%d\n", n(), J(), t(), U(), beta(), L());
}

void Par::print_label(FILE *file) const {
	fprintf(file, "n   J   t   U   beta   L\n");
}

// Var

void Var::print(FILE *file) const{
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
}

void Var::print_verboten(FILE *file) const{
	fprintf(file, "VARIABLES:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
	fprintf(file, "\n");
}

void Var::print_label(FILE *file) const {
	fprintf(file, "    barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d");
}

// Res

void Res::set(gsl_multiroot_fsolver *s) {
	for(int i=0; i<n; i++) {
		(*V)[i] = gsl_vector_get (s->x, i);
	}
	this->glob = case_2D2t_SGA_tPrim0::glob;
}

void Res::print(FILE *file) const { V->print(file); fprintf(file, "%21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f\n", glob.Energy, glob.gt, glob.gs, glob.gBDelta, glob.gm, glob.lbarchiAB, glob.lsimchiAB, glob.lDeltaS, glob.lDeltaT, glob.lsimDeltaB, glob.ln, glob.lmAF); }
void Res::print_verboten(FILE *file) const { fprintf(file, "Result: "); V->print_verboten(file); fprintf(file, "%21.15f %21.15f %21.15f %21.15f %21.15f\n", glob.Energy, glob.gt, glob.gs, glob.gBDelta, glob.gm); }
void Res::print_label(FILE *file) const { V->print_label(file); fprintf(file, "                      E                      gt                    gs             gBDelta               gm                 lbarchiAB                 lsimchiAB             lDeltaS             lDeltaT              lsimDeltaB            ln                  lmAF\n");}

void example() {

	printf("\nEXAMPLE 2D2t_tPrim0\n");

    Par par(N_PAR);
    par.set_n(0.9998);
    par.set_J(1./3.);
    par.set_t(-1.);
    par.set_U(11.2);
    par.set_beta(500.);
    par.set_L(1024);
    par.print_verboten();

    double set_barchiAB(3.362494837648e-01);
    double set_barDelta(3.350551731900e-01);
    double set_mu(3.550356000366e+00);
    double set_mAF(7.339835921625e-02);
    double set_d(6.726979905711e-05);
    double zero = 1e-6;

    Var initV(N_VAR);
    initV.set_barchiAB(0.5*set_barchiAB);
    initV.set_simchiAB(zero);
    initV.set_DeltaS(0.5*set_barDelta);
    initV.set_DeltaT(1e-12);
    initV.set_simDeltaB(zero);
    initV.set_mu(set_mu);
    initV.set_mAF(2.*set_mAF);
    initV.set_d(sqrt(set_d));

    // Z 1D1t_v2: U=10, n=0.97
	//sdelta          barchiAB             barDelta                mu                     mAF                    d                      E
	//                 3.119786e-01         -3.807814e-23          4.999908e+00         -1.624410e-17          4.852439e-01             4.848583	// 0, silne d
	//                -3.918124e-01          2.251022e-01          3.970843e+00         -2.576902e-11          2.888477e-02            -0.368634	// D
	//                -3.859593e-01          1.206862e-01          3.804065e+00         -1.879365e-01          3.084875e-02            -0.363913	// D, m
	//                -4.051486e-01          1.318796e-15          4.143907e+00          7.110104e-13          3.659600e-02            -0.353778	// 0
    // E:
    //11.200000 0.999800    3.362494837648e-01    3.350551731900e-01    3.550356000366e+00    7.339835921625e-02    6.726979905711e-05   -0.230396672929960     0.001202721086514     3.830458471110253     0.001208435395048     1.957155709469804     0.162201431523197     0.160426865804345    -3.550517931441191     0.110931280330479


    printf("Starting point: "); initV.print(); printf("\n\n");
    Res res(N_VAR); res.print_label();
    RootGSL root(par, function, res);
    if(root.Solve(initV) == GSL_SUCCESS) { root.print(); }

    printf("Koniec\n");
}

void example_1D_sdelta_run(Par &par, Var &initV) {
	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	initV.shake();
	if(root.Solve(initV) != GSL_SUCCESS) { printf("#D# First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#n           "); res.print_label();
	fflush(stdout);

	LOOPS = 20;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_n(0.998 + (double)i*0.0001);
		printf("%f ", par.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 348;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_n(0.998 - (double)i*0.001);
		printf("%f ", par.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
}

void example_1D_U_run(Par &par, Var &initV) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	int LOOPS;

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      "); res.print_label();
	fflush(stdout);

	initV.shake();
	if(root->Solve(initV) != GSL_SUCCESS) { printf("# First shot is not good.\n"); return; }
	Var old(*res.get_V());

	LOOPS = 100;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_U(10 - (double)i*0.1);
		printf("%f ", par.U());
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); continue; }
		root->print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 100;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_U(10 + (double)i*0.1);
		printf("%f ", par.U());
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); continue; }
		root->print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
	delete root;
}

void example_1D() {
    Par par(N_PAR);
    par.set_n(0.998);
    par.set_J(1./3.);
    par.set_t(-1.);
    par.set_U(12);
    par.set_beta(500.);
    par.set_L(1024);
    par.print_verboten();

	// I'm choosing start point
    Var initV_D(N_VAR);			// DS
    Var initV_E(N_VAR);			// m+DS+DT
    double ZERO = 1e-7;

//# mros-v42 - rozwazam model t-J-U z podwojna delta
//PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, U=11.200000, beta=500.000000, L=512
//sdelta          barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm
// hd      -1.747508266883526e-01 -2.027030520982555e-37 -1.240234181014150e-34 -9.718816970402378e-37 1.924280269183598e-40 5.599989556183054e+00 2.549859240836161e-17 7.059361913818627e-01     5.576444611323431    -0.002769198086339     0.000001914479987    -2.355321871631269    -0.831758073656629
// m+simD   2.024934059378212e-01 1.520174420318121e-29 2.437358465712652e-17 -3.704818546214106e-16 2.914037576947778e-02 5.324480009689015e+00 1.065093697377938e-03 1.425447264924440e-01      -0.198455179059096     0.170241383158736     3.650974765559812     0.822040769196572     1.421842230533131
// 0        2.026363818427083e-01 9.450840737303611e-28 -5.812953731560821e-17 -3.483503044169554e-17 -1.345335597537803e-17 5.326600208153319e+00 -2.645327197173785e-14 1.431820243635518e-01   -0.198343147089848     0.171570484014186     3.648200054185970     0.824185153602261     1.419294611259573
// hm+DS+DT 1.774033265920935e-01 1.510724295532038e-27 7.008119005952000e-02 2.511414612676356e-03 -7.644612066023643e-21 4.958059680002164e+00 5.544453294182708e-01 1.039146064012077e-01      -0.228849603677021     0.109295749934552     2.221902231754912     0.737471491338930     1.205713066275905
// DS       1.746930131714144e-01 -6.513397308362710e-30 1.634530545313107e-01 -3.917839824004685e-16 -2.694667188070913e-25 3.782641562157520e+00 -7.912028171777763e-14 4.061185417191907e-02   -0.244334759941973     0.026611929537691     3.941950059102869     0.311675662457952     1.829581856622030
// m+DS+DT  1.745312343776800e-01 -1.516038062635112e-25 1.632058553614935e-01 6.710615038101641e-05 -3.511378907752626e-23 3.771882611815728e+00 5.365718830046698e-02 4.049102383640613e-02     -0.244340250625107     0.026520366173390     3.919503626982180     0.310875778710769     1.824811298636322

    initV_D.set_barchiAB(1.746930131714144e-01);
    initV_D.set_simchiAB(ZERO);
    initV_D.set_DeltaS(1.634530545313107e-01);
    initV_D.set_DeltaT(ZERO);
    initV_D.set_simDeltaB(ZERO);
    initV_D.set_mu(3.782641562157520e+00);
    initV_D.set_mAF(ZERO);
    initV_D.set_d(4.061185417191907e-02);

    initV_E.set_barchiAB(1.745312343776800e-01);
    initV_E.set_simchiAB(ZERO);
    initV_E.set_DeltaS(1.632058553614935e-01);
    initV_E.set_DeltaT(6.710615038101641e-05);
    initV_E.set_simDeltaB(ZERO);
    initV_E.set_mu(3.771882611815728e+00);
    initV_E.set_mAF(5.365718830046698e-02);
    initV_E.set_d(4.049102383640613e-02);

//    example_1D_sdelta_run(par, initV_D);
    example_1D_sdelta_run(par, initV_E);
}

void example_2D_U_run(Par &par, Var &initV, int s1, int s2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	int LOOPS_n, LOOPS_U;
	Var initV_cp(N_VAR); initV_cp = initV;

	double start_U = 11.200000;
	double start_n = 0.996;

	par.set_U(start_U);
	par.set_n(start_n);

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta  U       "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	if(s2>0) { LOOPS_U = 188; } else { LOOPS_U = 112; }
	if(s1>0) { LOOPS_n = 40; } else { LOOPS_n = 60; }

	for(int i=0; i<=LOOPS_n; i++) {
		if(s1>0) { par.set_n(start_n + (double)i*0.0001); } else { par.set_n(start_n - (double)i*0.0001); }
		par.set_U(start_U);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS_U; j++) {
			if(s2>0) { par.set_U(start_U + (double)j*0.1); } else { par.set_U(start_U - (double)j*0.1); }
			printf("%f %f", par.n(), par.U());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}

		printf("\n");
		fflush(stdout);
	}
	printf("\n#END\n");
	fflush(stdout);
	initV = initV_cp;
}

void example_2D_n_run(Par &par, Var &initV, int s1, int s2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	int LOOPS_n, LOOPS_U;
	Var initV_cp(N_VAR); initV_cp = initV;

	double start_U = 11.200000;
	double start_n = 0.996;

	par.set_U(start_U);
	par.set_n(start_n);

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#U       sdelta        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	if(s2>0) { LOOPS_n = 0; } else { LOOPS_n = 198; }
	if(s1>0) { LOOPS_U = 88; } else { LOOPS_U = 112; }

	for(int i=0; i<=LOOPS_U; i++) {
		if(s1>0) { par.set_U(start_U + (double)i*0.1); } else { par.set_U(start_U - (double)i*0.1); }
		par.set_n(start_n);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS_n; j++) {
			if(s2>0) { par.set_n(start_n + (double)j*0.002); } else { par.set_n(start_n - (double)j*0.002); }
			printf("%f %f", par.U(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}

		printf("\n");
		fflush(stdout);
	}
	printf("\n#END\n");
	fflush(stdout);
	initV = initV_cp;
}

void example_2D_Jn_run(Par &par, Var &initV, int s1, int s2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	int LOOPS_n, LOOPS_J;
	Var initV_cp(N_VAR); initV_cp = initV;

	double start_J = 0.33;
	double start_n = 0.999;

	par.set_J(start_J);
	par.set_n(start_n);

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#J       sdelta        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	if(s1<0) { LOOPS_J = 33; } else { LOOPS_J = 67; }
	if(s2<0) { LOOPS_n = 60; } else { LOOPS_n = 10; }

	for(int i=0; i<=LOOPS_J; i++) {
		if(s1<0) { par.set_J(start_J - (double)i*0.01); } else { par.set_J(start_J + (double)i*0.01); }
		par.set_n(start_n);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS_n; j++) {
			if(s2<0) { par.set_n(start_n - (double)j*0.0001); } else { par.set_n(start_n + (double)j*0.0001); }
			printf("%f %f", par.J(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}

		printf("\n");
		fflush(stdout);
	}
	printf("\n#END\n");
	fflush(stdout);
	initV = initV_cp;
}

void example_2D() {
	Par par(N_PAR);
	par.set_n(0.999);
	par.set_J(1./3.);
	par.set_t(-1.);
	par.set_U(100);
	par.set_beta(1500);
	par.set_L(512);
	printf("#"); par.print_verboten(); fflush(stdout);

	// I'm choosing start point
    Var initV_D(N_VAR);			// DS
    Var initV_E(N_VAR);			// m+DS+DT
    double ZERO = 1e-7;

//# mros-v42 - rozwazam model t-J-U z podwojna delta
//PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, U=11.200000, beta=500.000000, L=512
//sdelta          barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm
// hd      -1.747508266883526e-01 -2.027030520982555e-37 -1.240234181014150e-34 -9.718816970402378e-37 1.924280269183598e-40 5.599989556183054e+00 2.549859240836161e-17 7.059361913818627e-01     5.576444611323431    -0.002769198086339     0.000001914479987    -2.355321871631269    -0.831758073656629
// m+simD   2.024934059378212e-01 1.520174420318121e-29 2.437358465712652e-17 -3.704818546214106e-16 2.914037576947778e-02 5.324480009689015e+00 1.065093697377938e-03 1.425447264924440e-01      -0.198455179059096     0.170241383158736     3.650974765559812     0.822040769196572     1.421842230533131
// 0        2.026363818427083e-01 9.450840737303611e-28 -5.812953731560821e-17 -3.483503044169554e-17 -1.345335597537803e-17 5.326600208153319e+00 -2.645327197173785e-14 1.431820243635518e-01   -0.198343147089848     0.171570484014186     3.648200054185970     0.824185153602261     1.419294611259573
// hm+DS+DT 1.774033265920935e-01 1.510724295532038e-27 7.008119005952000e-02 2.511414612676356e-03 -7.644612066023643e-21 4.958059680002164e+00 5.544453294182708e-01 1.039146064012077e-01      -0.228849603677021     0.109295749934552     2.221902231754912     0.737471491338930     1.205713066275905
// DS       1.746930131714144e-01 -6.513397308362710e-30 1.634530545313107e-01 -3.917839824004685e-16 -2.694667188070913e-25 3.782641562157520e+00 -7.912028171777763e-14 4.061185417191907e-02   -0.244334759941973     0.026611929537691     3.941950059102869     0.311675662457952     1.829581856622030
// m+DS+DT  1.745312343776800e-01 -1.516038062635112e-25 1.632058553614935e-01 6.710615038101641e-05 -3.511378907752626e-23 3.771882611815728e+00 5.365718830046698e-02 4.049102383640613e-02     -0.244340250625107     0.026520366173390     3.919503626982180     0.310875778710769     1.824811298636322

    initV_D.set_barchiAB(1.746930131714144e-01);
    initV_D.set_simchiAB(ZERO);
    initV_D.set_DeltaS(1.634530545313107e-01);
    initV_D.set_DeltaT(ZERO);
    initV_D.set_simDeltaB(ZERO);
    initV_D.set_mu(3.782641562157520e+00);
    initV_D.set_mAF(ZERO);
    initV_D.set_d(4.061185417191907e-02);

    initV_E.set_barchiAB(1.745312343776800e-01);
    initV_E.set_simchiAB(ZERO);
    initV_E.set_DeltaS(1.632058553614935e-01);
    initV_E.set_DeltaT(6.710615038101641e-05);
    initV_E.set_simDeltaB(ZERO);
    initV_E.set_mu(3.771882611815728e+00);
    initV_E.set_mAF(5.365718830046698e-02);
    initV_E.set_d(4.049102383640613e-02);

//    	printf("# CASE D SGA_v2++\n"); example_2D_U_run(par, initV_D, 1, 1);
//    	printf("# CASE D SGA_v2+-\n"); example_2D_U_run(par, initV_D, 1, -1);
//    	printf("# CASE D SGA_v2-+\n"); example_2D_U_run(par, initV_D, -1, 1);
//    	printf("# CASE D SGA_v2--\n"); example_2D_U_run(par, initV_D, -1, -1);

//      	  printf("# CASE E SGA_v2++\n"); example_2D_U_run(par, initV_E, 1, 1);
//      	  printf("# CASE E SGA_v2+-\n"); example_2D_U_run(par, initV_E, 1, -1);
//      	  printf("# CASE E SGA_v2-+\n"); example_2D_U_run(par, initV_E, -1, 1);
//      	  printf("# CASE E SGA_v2--\n"); example_2D_U_run(par, initV_E, -1, -1);

//    	printf("# CASE D SGA_v2++\n"); example_2D_n_run(par, initV_D, 1, 1);
//    	printf("# CASE D SGA_v2+-\n"); example_2D_n_run(par, initV_D, 1, -1);
//    	printf("# CASE D SGA_v2-+\n"); example_2D_n_run(par, initV_D, -1, 1);
//    	printf("# CASE D SGA_v2--\n"); example_2D_n_run(par, initV_D, -1, -1);


    //	printf("# CASE E SGA_v2+-\n"); example_2D_n_run(par, initV_E, 1, -1);
    	//printf("# CASE E SGA_v2-+\n"); example_2D_n_run(par, initV_E, -1, 1);
		//printf("# CASE E SGA_v2--\n"); example_2D_n_run(par, initV_E, -1, -1);
    	//printf("# CASE E SGA_v2++\n"); example_2D_n_run(par, initV_E, 1, 1);

    // #sdelta  U           barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm
    //0.996000 30.000000 1.712339149661826e-01 -7.116046177010064e-25 1.667588557806268e-01 8.099181961153371e-05 2.475126730413155e-24 2.543659754840674e+00 6.868129192639835e-02 6.918257610033627e-03   -0.239724543090433     0.009881491279315     3.930261184553631     0.009941314187715     1.982488634154969
    Var initV_U30_J13_delta996(N_VAR);	// m+DS+DT
    initV_U30_J13_delta996.set_barchiAB(1.712339149661826e-01);
    initV_U30_J13_delta996.set_simchiAB(ZERO);
    initV_U30_J13_delta996.set_DeltaS(1.667588557806268e-01);
    initV_U30_J13_delta996.set_DeltaT(8.099181961153371e-05);
    initV_U30_J13_delta996.set_simDeltaB(ZERO);
    initV_U30_J13_delta996.set_mu(2.543659754840674e+00);
    initV_U30_J13_delta996.set_mAF(6.868129192639835e-02);
    initV_U30_J13_delta996.set_d(6.918257610033627e-03);

    Var initV_U100_J13_delta999(N_VAR);	// m+DS+DT
    initV_U100_J13_delta999.set_barchiAB(1.686360066445106e-01);
    initV_U100_J13_delta999.set_simchiAB(ZERO);
    initV_U100_J13_delta999.set_DeltaS(1.676171218402892e-01);
    initV_U100_J13_delta999.set_DeltaT(3.872342472544926e-05);
    initV_U100_J13_delta999.set_simDeltaB(ZERO);
    initV_U100_J13_delta999.set_mu(2.266511470421953e+00);
    initV_U100_J13_delta999.set_mAF(1.319374026200395e-01);
    initV_U100_J13_delta999.set_d(9.082992468738691e-04);

    //printf("# CASE EJ SGA_v2++\n"); example_2D_Jn_run(par, initV_U30_J13_delta996, 1, 1);
    //printf("# CASE EJ SGA_v2+-\n"); example_2D_Jn_run(par, initV_U30_J13_delta996, 1, -1);
    printf("# CASE EJ SGA_v2-+\n"); example_2D_Jn_run(par, initV_U100_J13_delta999, -1, 1);
    printf("# CASE EJ SGA_v2--\n"); example_2D_Jn_run(par, initV_U100_J13_delta999, -1, -1);
}

void example_inf() {
    Par par(N_PAR);
    par.set_n(0.996);
    par.set_J(1./3.);
    par.set_t(-1.);
    par.set_U(8);
    par.set_beta(500.);
    par.set_L(256);
    par.print_verboten();

	// I'm choosing start point
	Var initV(N_VAR);

	FILE * pFile;
	pFile = fopen(OUTPUT_FILE, "w");

	// INFINITE
	Res res(N_VAR);
	RootGSL root(par, function, res);
	printf("sdelta      "); res.print_label();
	while(1) {
	    initV.set_barchiAB(frand());
	    initV.set_simchiAB(frand());
	    initV.set_DeltaS(frand());
	    initV.set_DeltaT(frand());
	    initV.set_simDeltaB(frand());
	    initV.set_mu(16.*frand()-8.);
	    initV.set_mAF(frand());
	    initV.set_d(frand());

		if(root.Solve(initV) == GSL_SUCCESS) {
			printf("       "); root.print(); root.print(pFile); fflush(pFile); fflush(stdout);
		}
	}
	fclose (pFile);
}

void test_energy() {};

void run_1D_n() {
    Par par(N_PAR);
    par.set_n(0.996);
    par.set_J(1./3.);
    par.set_t(-1.);
    par.set_U(12);
    par.set_beta(20.);
    par.set_L(512);
	printf("#"); par.print_label();
	printf("#"); par.print_verboten();

	//PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, U=11.200000, beta=500.000000, L=512
	//sdelta          barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm

	Var initV_A(N_VAR); // 0
	Var initV_B(N_VAR); // m+simD
	Var initV_C(N_VAR); // DS
	Var initV_D(N_VAR);	// hm+DS+DT
    Var initV_E(N_VAR);	// m+DS+DT
    Var initV_E2(N_VAR);	// m+DS+DT
    double ZERO = 1e-7;

    //12.000000  2.026184381598744e-01 -1.186418777640084e-27 -3.429444614835553e-24 1.526293500665240e-23 -1.057524786505130e-23 5.298201527290678e+00 2.190787957523977e-13 8.516497617012024e-02   -0.187583720791968     0.071797097099770     3.853444048374414     0.071797097099770     1.963019115641622
    initV_A.set_barchiAB(2.026184381598744e-01);
    initV_A.set_simchiAB(ZERO);
    initV_A.set_DeltaS(ZERO);
    initV_A.set_DeltaT(ZERO);
    initV_A.set_simDeltaB(ZERO);
    initV_A.set_mu(5.298201527290678e+00);
    initV_A.set_mAF(ZERO);
    initV_A.set_d(8.516497617012024e-02);

    //12.000000  2.021237720422460e-01 -5.244786255648164e-29 -1.738944248870892e-26 -1.236739598082204e-26 4.720333519214653e-02 5.267184394745283e+00 9.338566671508935e-04 8.289336108866857e-02   -0.187965111416371     0.068815038617225     3.859435022562903     0.068815040074678     1.964544482205201
    initV_B.set_barchiAB(2.021237720422460e-01);
    initV_B.set_simchiAB(ZERO);
    initV_B.set_DeltaS(ZERO);
    initV_B.set_DeltaT(ZERO);
    initV_B.set_simDeltaB(4.720333519214653e-02);
    initV_B.set_mu(5.267184394745283e+00);
    initV_B.set_mAF(9.338566671508935e-04);
    initV_B.set_d(8.289336108866857e-02);

    //12.000000  1.736133820319831e-01 2.629411885476526e-34 1.648124407812186e-01 2.180115868042694e-13 1.508996038778887e-32 3.503145875806553e+00 2.249459583713295e-09 3.132854735959725e-02   -0.243400155333143     0.020647011628417     3.952565062361090     0.020647011628417     1.988105898175721
    initV_C.set_barchiAB(1.736133820319831e-01);
    initV_C.set_simchiAB(ZERO);
    initV_C.set_DeltaS(1.648124407812186e-01);
    initV_C.set_DeltaT(ZERO);
    initV_C.set_simDeltaB(ZERO);
    initV_C.set_mu(3.503145875806553e+00);
    initV_C.set_mAF(ZERO);
    initV_C.set_d(3.132854735959725e-02);

    //12.000000  1.698382752577751e-01 3.924277188550582e-28 5.087494566443462e-02 3.915321475979912e-03 2.031881214632224e-26 4.924395250998753e+00 6.357958239182453e-01 8.209765805706415e-02   -0.222097059328138     0.080590383799718     1.958268476595549     0.081372214625683     1.399381462145168
    initV_D.set_barchiAB(1.698382752577751e-01);
    initV_D.set_simchiAB(ZERO);
    initV_D.set_DeltaS(5.087494566443462e-02);
    initV_D.set_DeltaT(3.915321475979912e-03);
    initV_D.set_simDeltaB(ZERO);
    initV_D.set_mu(4.924395250998753e+00);
    initV_D.set_mAF(6.357958239182453e-01);
    initV_D.set_d(8.209765805706415e-02);

    //12.000000  1.733764351054193e-01 7.651187477648849e-28 1.644567593975851e-01 7.941687221654041e-05 2.015988158887319e-27 3.488839270407357e+00 6.465659750128631e-02 3.125578512762359e-02   -0.243411142424865     0.020592193200615     3.919795243819099     0.020617325465607     1.979847277902793
    initV_E.set_barchiAB(1.733764351054193e-01);
    initV_E.set_simchiAB(ZERO);
    initV_E.set_DeltaS(1.644567593975851e-01);
    initV_E.set_DeltaT(7.941687221654041e-05);
    initV_E.set_simDeltaB(ZERO);
    initV_E.set_mu(3.488839270407357e+00);
    initV_E.set_mAF(6.465659750128631e-02);
    initV_E.set_d(3.125578512762359e-02);

    // U = 5.0

    //5.000000  2.022852660351944e-01 8.314687428172392e-29 4.114238988963787e-02 1.004077165481028e-13 -7.685333185231544e-27 2.475473561708106e+00 2.913553012463496e-12 3.769465553943369e-01   -0.703722356482991     0.820498368292789     2.026830449489831     0.820498368292788     1.423667956192676     0.922997852463609     0.000000000000000     0.020847162148275    -0.000000000000017     0.000000000000000    -4.960045018989471     0.000000000001923
    initV_E2.set_barchiAB(2.022852660351944e-01);
    initV_E2.set_simchiAB(ZERO);
    initV_E2.set_DeltaS(4.114238988963787e-02);
    initV_E2.set_DeltaT(ZERO);
    initV_E2.set_simDeltaB(ZERO);
    initV_E2.set_mu(2.475473561708106e+00);
    initV_E2.set_mAF(ZERO);
    initV_E2.set_d(3.769465553943369e-01);

    // (0.004,12) G "rowniez D m and BDelta"
    // #sdelta  U           barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm
    //0.996000 30.000000 1.712339149661826e-01 -7.116046177010064e-25 1.667588557806268e-01 8.099181961153371e-05 2.475126730413155e-24 2.543659754840674e+00 6.868129192639835e-02 6.918257610033627e-03   -0.239724543090433     0.009881491279315     3.930261184553631     0.009941314187715     1.982488634154969
    Var initV_U30(N_VAR);	// m+DS+DT
    initV_U30.set_barchiAB(1.712339149661826e-01);
    initV_U30.set_simchiAB(ZERO);
    initV_U30.set_DeltaS(1.667588557806268e-01);
    initV_U30.set_DeltaT(8.099181961153371e-05);
    initV_U30.set_simDeltaB(ZERO);
    initV_U30.set_mu(2.543659754840674e+00);
    initV_U30.set_mAF(6.868129192639835e-02);
    initV_U30.set_d(6.918257610033627e-03);

    //12.000000  1.733802388850798e-01 -6.891210342912008e-26 1.643694744620512e-01 8.105571459244505e-05 -2.390068851350686e-25 3.489032035238794e+00 6.592870961037603e-02 3.126274297541825e-02   -0.243337347752307     0.020595761863424     3.918491801453419     0.020621887467297     1.979518073030256     0.190443023014729    -0.000000000000000     0.161020109522189    -0.000026468012758     0.000000000000000    -6.984734776779710     0.052223578354371
    Var initV_U12(N_VAR);	// m+DS+DT
    initV_U12.set_barchiAB(1.733764351054193e-01);
    initV_U12.set_simchiAB(ZERO);
    initV_U12.set_DeltaS(1.644567593975851e-01);
    initV_U12.set_DeltaT(8.105571459244505e-05);
    initV_U12.set_simDeltaB(ZERO);
    initV_U12.set_mu(3.488839270407357e+00);
    initV_U12.set_mAF(6.465659750128631e-02);
    initV_U12.set_d(3.125578512762359e-02);

    //8.000000  2.005663179081676e-01 -2.904629914076084e-28 7.776740184343016e-02 4.050238821696296e-14 2.367150445026661e-26 3.937985231588939e+00 1.070519194980411e-10 2.831618246729351e-01   -0.369202945130017     0.549313915936231     2.793257946407135     0.549313915936231     1.671304265059817     0.689372281255883    -0.000000000000000     0.054306103292650    -0.000000000000009    -0.000000000000000    -7.887674981155220     0.000000000088600
    Var initV_U8(N_VAR);	// m+DS+DT
    initV_U8.set_barchiAB(2.005663179081676e-01);
    initV_U8.set_simchiAB(ZERO);
    initV_U8.set_DeltaS(7.776740184343016e-02);
    initV_U8.set_DeltaT(ZERO);
    initV_U8.set_simDeltaB(ZERO);
    initV_U8.set_mu(3.937985231588939e+00);
    initV_U8.set_mAF(ZERO);
    initV_U8.set_d(2.831618246729351e-01);


    //A 0
    //B m+simD
    //C DS
    //D hm+DS+DT
    //E m+DS+DT
	// I'm choosing start point
	Var initV = initV_E;

	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      "); res.print_label();
	fflush(stdout);

	int MNOZNIK = 1;

	LOOPS = 40*MNOZNIK;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_n(0.996 + ((double)i)*0.0001/MNOZNIK);
		printf("%f ", 1.-par.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 60*MNOZNIK;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_n(0.996 - ((double)i)*0.0001/MNOZNIK);
		printf("%f ", 1.-par.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET\n");

	LOOPS = 590*MNOZNIK;
	for(int i=0; i<=LOOPS; i++) {
		par.set_n(0.99 - ((double)i)*0.001/MNOZNIK);
		printf("%f ", 1.-par.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
};

void run_1D_U() {
    Par par(N_PAR);
    par.set_n(0.95);
    par.set_J(1./3.);
    par.set_t(-1.);
    par.set_U(30);
    par.set_beta(1500.);
    par.set_L(512);
	printf("#"); par.print_label();
	printf("#"); par.print_verboten();

	//PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, U=11.200000, beta=500.000000, L=512
	//sdelta          barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm
	// hd      -1.747508266883526e-01 -2.027030520982555e-37 -1.240234181014150e-34 -9.718816970402378e-37 1.924280269183598e-40 5.599989556183054e+00 2.549859240836161e-17 7.059361913818627e-01     5.576444611323431    -0.002769198086339     0.000001914479987    -2.355321871631269    -0.831758073656629
	// m+simD   2.024934059378212e-01 1.520174420318121e-29 2.437358465712652e-17 -3.704818546214106e-16 2.914037576947778e-02 5.324480009689015e+00 1.065093697377938e-03 1.425447264924440e-01      -0.198455179059096     0.170241383158736     3.650974765559812     0.822040769196572     1.421842230533131
	// 0        2.026363818427083e-01 9.450840737303611e-28 -5.812953731560821e-17 -3.483503044169554e-17 -1.345335597537803e-17 5.326600208153319e+00 -2.645327197173785e-14 1.431820243635518e-01   -0.198343147089848     0.171570484014186     3.648200054185970     0.824185153602261     1.419294611259573
	// hm+DS+DT 1.774033265920935e-01 1.510724295532038e-27 7.008119005952000e-02 2.511414612676356e-03 -7.644612066023643e-21 4.958059680002164e+00 5.544453294182708e-01 1.039146064012077e-01      -0.228849603677021     0.109295749934552     2.221902231754912     0.737471491338930     1.205713066275905
	// DS       1.746930131714144e-01 -6.513397308362710e-30 1.634530545313107e-01 -3.917839824004685e-16 -2.694667188070913e-25 3.782641562157520e+00 -7.912028171777763e-14 4.061185417191907e-02   -0.244334759941973     0.026611929537691     3.941950059102869     0.311675662457952     1.829581856622030
	// m+DS+DT  1.745312343776800e-01 -1.516038062635112e-25 1.632058553614935e-01 6.710615038101641e-05 -3.511378907752626e-23 3.771882611815728e+00 5.365718830046698e-02 4.049102383640613e-02     -0.244340250625107     0.026520366173390     3.919503626982180     0.310875778710769     1.824811298636322

	Var initV_A(N_VAR); // 0
	Var initV_B(N_VAR); // m+simD
	Var initV_C(N_VAR); // DS
	Var initV_D(N_VAR);	// hm+DS+DT
    Var initV_E(N_VAR);	// m+DS+DT
    double ZERO = 1e-7;

    initV_A.set_barchiAB(2.026363818427083e-01);
    initV_A.set_simchiAB(ZERO);
    initV_A.set_DeltaS(ZERO);
    initV_A.set_DeltaT(ZERO);
    initV_A.set_simDeltaB(ZERO);
    initV_A.set_mu(5.326600208153319e+00);
    initV_A.set_mAF(ZERO);
    initV_A.set_d(1.431820243635518e-01);

    initV_B.set_barchiAB(2.024934059378212e-01);
    initV_B.set_simchiAB(ZERO);
    initV_B.set_DeltaS(ZERO);
    initV_B.set_DeltaT(ZERO);
    initV_B.set_simDeltaB(2.914037576947778e-02);
    initV_B.set_mu(5.324480009689015e+00);
    initV_B.set_mAF(1.065093697377938e-03);
    initV_B.set_d(1.425447264924440e-01);

    initV_C.set_barchiAB(1.746930131714144e-01);
    initV_C.set_simchiAB(ZERO);
    initV_C.set_DeltaS(1.634530545313107e-01);
    initV_C.set_DeltaT(ZERO);
    initV_C.set_simDeltaB(ZERO);
    initV_C.set_mu(3.782641562157520e+00);
    initV_C.set_mAF(ZERO);
    initV_C.set_d(4.061185417191907e-02);

    initV_D.set_barchiAB(1.774033265920935e-01);
    initV_D.set_simchiAB(ZERO);
    initV_D.set_DeltaS(7.008119005952000e-02);
    initV_D.set_DeltaT(2.511414612676356e-03);
    initV_D.set_simDeltaB(ZERO);
    initV_D.set_mu(4.958059680002164e+00);
    initV_D.set_mAF(5.544453294182708e-01);
    initV_D.set_d(1.039146064012077e-01);

    // m+DS+DT  1.745312343776800e-01 -1.516038062635112e-25 1.632058553614935e-01 6.710615038101641e-05 -3.511378907752626e-23 3.771882611815728e+00 5.365718830046698e-02 4.049102383640613e-02     -0.244340250625107     0.026520366173390     3.919503626982180     0.310875778710769     1.824811298636322
    initV_E.set_barchiAB(1.745312343776800e-01);
    initV_E.set_simchiAB(ZERO);
    initV_E.set_DeltaS(1.632058553614935e-01);
    initV_E.set_DeltaT(6.710615038101641e-05);
    initV_E.set_simDeltaB(ZERO);
    initV_E.set_mu(3.771882611815728e+00);
    initV_E.set_mAF(5.365718830046698e-02);
    initV_E.set_d(4.049102383640613e-02);

    // U=12
	//#sdelta          barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm                 lbarchiAB                 lsimchiAB             lDeltaS             lDeltaT              lsimDeltaB            ln                  lmAF
    //0.000000  1.674082457211893e-01 -2.093294063087793e-28 4.469384019075074e-02 7.929432187697337e-14 -1.733965103873139e-26 5.999999999998983e+00 6.573393113098631e-01 8.596995357628882e-02   -0.200409018728056     0.071625608076832     1.893133858068533     0.071625608076832     1.375912009566213     0.150857162600492    -0.000000000000000     0.021152855528054    -0.000000000000013     0.000000000000000   -11.999999999999901     0.388287452468716
    Var initV_E3(N_VAR);	// m+DS+DT
    initV_E3.set_barchiAB(1.674082457211893e-01);
    initV_E3.set_simchiAB(ZERO);
    initV_E3.set_DeltaS(4.469384019075074e-02);
    initV_E3.set_DeltaT(7.929432187697337e-14);
    initV_E3.set_simDeltaB(ZERO);
    initV_E3.set_mu(5.999999999998983e+00);
    initV_E3.set_mAF(6.573393113098631e-01);
    initV_E3.set_d(8.596995357628882e-02);

    //0.001000  1.693211838431733e-01 -9.010683092846660e-27 1.669113860227631e-01 3.921232154000801e-05 1.795614877775459e-27 3.371546654213554e+00 1.319009092615469e-01 1.496850257060662e-02   -0.233035552832709     0.004971366186468     3.853190565855477     0.004999079730841     1.962954550124754     0.168078063232467    -0.000000000000000     0.160785344489193    -0.000012591045619    -0.000000000000000    -6.744734991792987     0.100877447601264
    Var initV_E4(N_VAR);	// m+DS+DT
    initV_E4.set_barchiAB(1.693211838431733e-01);
    initV_E4.set_simchiAB(ZERO);
    initV_E4.set_DeltaS(1.669113860227631e-01);
    initV_E4.set_DeltaT(3.921232154000801e-05);
    initV_E4.set_simDeltaB(ZERO);
    initV_E4.set_mu(3.371546654213554e+00);
    initV_E4.set_mAF(1.319009092615469e-01);
    initV_E4.set_d(1.496850257060662e-02);

    // #sdelta  U           barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm
    //0.996000 30.000000 1.712339149661826e-01 -7.116046177010064e-25 1.667588557806268e-01 8.099181961153371e-05 2.475126730413155e-24 2.543659754840674e+00 6.868129192639835e-02 6.918257610033627e-03   -0.239724543090433     0.009881491279315     3.930261184553631     0.009941314187715     1.982488634154969
    Var initV_n996(N_VAR);	// m+DS+DT
    initV_n996.set_barchiAB(1.712339149661826e-01);
    initV_n996.set_simchiAB(ZERO);
    initV_n996.set_DeltaS(1.667588557806268e-01);
    initV_n996.set_DeltaT(8.099181961153371e-05);
    initV_n996.set_simDeltaB(ZERO);
    initV_n996.set_mu(2.543659754840674e+00);
    initV_n996.set_mAF(6.868129192639835e-02);
    initV_n996.set_d(6.918257610033627e-03);

    //0.000100  1.678348042450931e-01 -2.927395360969979e-28 1.677167865222651e-01 4.376891846493037e-06 -1.221617979160505e-27 2.464270634969308e+00 1.491573220141453e-01 1.088903265198284e-03   -0.229934330728995     0.000245072324020     3.827001361851074     0.000252233993712     1.956272312805933     0.160821078427015    -0.000000000000000     0.160463092606499    -0.000001395864255     0.000000000000000    -4.928703683852456     0.112372751624520
    Var initV_n9999(N_VAR);	// m+DS+DT
    initV_n9999.set_barchiAB(1.678348042450931e-01);
    initV_n9999.set_simchiAB(ZERO);
    initV_n9999.set_DeltaS(1.677167865222651e-01);
    initV_n9999.set_DeltaT(4.376891846493037e-06);
    initV_n9999.set_simDeltaB(ZERO);
    initV_n9999.set_mu(2.464270634969308e+00);
    initV_n9999.set_mAF(1.491573220141453e-01);
    initV_n9999.set_d(1.088903265198284e-03);

    //0.005000  1.719360391304059e-01 2.884269149770376e-27 1.664233461987727e-01 6.009383495201743e-05 8.069338879136127e-27 2.558033044487238e+00 4.084324677579575e-02 7.743651247825312e-03   -0.242260650547522     0.012366483821177     3.946166323212613     0.012392829750394     1.986495991240005     0.181988535661920     0.000000000000000     0.164183551041488    -0.000019761688977    -0.000000000000000    -5.124035767127451     0.031596149145143
    Var initV_n995(N_VAR);	// m+DS+DT
    initV_n995.set_barchiAB(1.719360391304059e-01);
    initV_n995.set_simchiAB(ZERO);
    initV_n995.set_DeltaS(1.664233461987727e-01);
    initV_n995.set_DeltaT(6.009383495201743e-05);
    initV_n995.set_simDeltaB(ZERO);
    initV_n995.set_mu(2.558033044487238e+00);
    initV_n995.set_mAF(4.084324677579575e-02);
    initV_n995.set_d(7.743651247825312e-03);

    //0.010000  1.743690276673767e-01 2.423962625931899e-28 1.638240257680628e-01 3.056249295120055e-14 9.615692361877309e-26 2.578367387013500e+00 4.016959150957403e-15 1.104740366571901e-02   -0.255116725490952     0.024681067531451     3.919250854474984     0.024681067531452     1.979709790468033     0.195530057701286     0.000000000000000     0.160516863243753    -0.000000000000010    -0.000000000000000    -5.172753486009316     0.000000000000003
    Var initV_n99(N_VAR);	// m+DS+DT
    initV_n99.set_barchiAB(1.743690276673767e-01);
    initV_n99.set_simchiAB(ZERO);
    initV_n99.set_DeltaS(1.638240257680628e-01);
    initV_n99.set_DeltaT(ZERO);
    initV_n99.set_simDeltaB(ZERO);
    initV_n99.set_mu(2.578367387013500e+00);
    initV_n99.set_mAF(ZERO);
    initV_n99.set_d(1.104740366571901e-02);

    //0.030000  1.813520231990706e-01 1.409966887506323e-26 1.535516973330122e-01 -2.461608352911993e-14 6.144942323143734e-26 2.587373612445779e+00 2.021733708709011e-14 1.950807742883326e-02   -0.306882724994984     0.072878021070937     3.764468942033635     0.072878021070937     1.940223941207209     0.243551535797903     0.000000000000000     0.144510148901668     0.000000000000008    -0.000000000000000    -5.228249747759156     0.000000000000011
    Var initV_n97(N_VAR);	// m+DS+DT
    initV_n97.set_barchiAB(1.813520231990706e-01);
    initV_n97.set_simchiAB(ZERO);
    initV_n97.set_DeltaS(1.535516973330122e-01);
    initV_n97.set_DeltaT(ZERO);
    initV_n97.set_simDeltaB(ZERO);
    initV_n97.set_mu(2.587373612445779e+00);
    initV_n97.set_mAF(ZERO);
    initV_n97.set_d(1.950807742883326e-02);

    //0.050000  1.859584750151545e-01 -1.759399038333364e-28 1.440054531663302e-01 2.446297580902336e-12 -2.862020734914765e-26 2.547992271171315e+00 -3.879912715565515e-14 2.522884331606858e-02   -0.358294445844847     0.119130166754930     3.618401155843670     0.119130166754930     1.902209545724043     0.287348256988370    -0.000000000000000     0.130267374546210    -0.000000000000738     0.000000000000000    -5.192521621059769    -0.000000000000014
    Var initV_n95(N_VAR);	// m+DS+DT
    initV_n95.set_barchiAB(1.859584750151545e-01);
    initV_n95.set_simchiAB(ZERO);
    initV_n95.set_DeltaS(1.440054531663302e-01);
    initV_n95.set_DeltaT(ZERO);
    initV_n95.set_simDeltaB(ZERO);
    initV_n95.set_mu(2.547992271171315e+00);
    initV_n95.set_mAF(ZERO);
    initV_n95.set_d(2.522884331606858e-02);

    //0.100000  1.922480249594002e-01 1.281934136381806e-27 1.216499929328693e-01 2.079431279827296e-15 -3.442946450801562e-26 2.355075889970666e+00 1.270330041753321e-15 3.458764721958481e-02   -0.481237536927904     0.225577291012653     3.288231916584824     0.225577291012653     1.813348261251772     0.383616313905627     0.000000000000000     0.100003347353545    -0.000000000000001     0.000000000000000    -4.930519076241818    -0.000000000000000
    Var initV_n9(N_VAR);	// m+DS+DT
    initV_n9.set_barchiAB(1.922480249594002e-01);
    initV_n9.set_simchiAB(ZERO);
    initV_n9.set_DeltaS(1.216499929328693e-01);
    initV_n9.set_DeltaT(ZERO);
    initV_n9.set_simDeltaB(ZERO);
    initV_n9.set_mu(2.355075889970666e+00);
    initV_n9.set_mAF(ZERO);
    initV_n9.set_d(3.458764721958481e-02);

	// I'm choosing start point
	Var initV = initV_n9;

	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#U           "); res.print_label();
	fflush(stdout);

	/*
	LOOPS = 0;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_U(30 + ((double)i)*0.1);
		printf("%f ", par.U());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");
	*/

	LOOPS = 300;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_U(30 - ((double)i)*0.1);
		printf("%f ", par.U());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
};

void test_L() {
	Par par(N_PAR);
	par.set_n(0.7);
	par.set_J(1./3.);
	par.set_t(-1.);
	par.set_U(5);
	par.set_beta(1500.);
	par.set_L(128);
	printf("#"); par.print_label();
	printf("#"); par.print_verboten();

	Var initV_A_n999(N_VAR); // 0
	Var initV_B_n999(N_VAR); // m+simD
	Var initV_C_n999(N_VAR); // DS
	Var initV_D_n999(N_VAR); // hm+DS+DT
    Var initV_E_n999(N_VAR); // m+DS+DT

	Var initV_A_n97(N_VAR); // 0
	Var initV_B_n97(N_VAR); // m+simD
	Var initV_C_n97(N_VAR); // DS
	Var initV_D_n97(N_VAR); // hm+DS+DT
    Var initV_E_n97(N_VAR); // m+DS+DT

    Var initV_E_n3(N_VAR); // m+DS+DT

	double ZERO = 1e-7;

	//#sdelta          barchiAB             simchiAB                 DeltaS             DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm                 lbarchiAB                 lsimchiAB             lDeltaS             lDeltaT              lsimDeltaB            ln                  lmAF

	// A
	// 0.001000  2.026327689968521e-01 2.516515649756660e-30 -8.014736079845175e-22 2.617949860434002e-22 3.519093047806009e-21 5.755090124958051e+00 -1.959379391656846e-11 7.770809162945296e-02   -0.171106337045195     0.051548809695281     3.896075073726727     0.051548809695282     1.973847783829018     0.248916929797497     0.000000000000000    -0.000000000000000    -0.000000000000000    -0.000000000000000   -11.510812906467720    -0.000000000025807
	// 0.030000  2.025644934728421e-01 1.361964670396322e-28 3.633137078445894e-24 4.595688949245551e-23 -2.188197897581085e-21 4.132763653906881e+00 7.300750618813632e-13 1.164429991364329e-01   -0.306859077264884     0.199548389452324     3.562515928417966     0.199548389452324     1.887462828353969     0.379958198084553     0.000000000000000     0.000000000000000    -0.000000000000000     0.000000000000000    -8.298225286901696     0.000000000000804

    initV_A_n999.set_barchiAB(2.026327689968521e-01);
    initV_A_n999.set_simchiAB(ZERO);
    initV_A_n999.set_DeltaS(ZERO);
    initV_A_n999.set_DeltaT(ZERO);
    initV_A_n999.set_simDeltaB(ZERO);
    initV_A_n999.set_mu(5.755090124958051e+00);
    initV_A_n999.set_mAF(ZERO);
    initV_A_n999.set_d(7.770809162945296e-02);

    initV_A_n97.set_barchiAB(2.025644934728421e-01);
    initV_A_n97.set_simchiAB(ZERO);
    initV_A_n97.set_DeltaS(ZERO);
    initV_A_n97.set_DeltaT(ZERO);
    initV_A_n97.set_simDeltaB(ZERO);
    initV_A_n97.set_mu(4.132763653906881e+00);
    initV_A_n97.set_mAF(ZERO);
    initV_A_n97.set_d(1.164429991364329e-01);

	// B
	// 0.001000  2.019508440767480e-01 -6.121550827870031e-28 2.201084548781910e-15 -5.079815375594085e-15 5.253586960997034e-02 5.721566801475077e+00 2.439528519551680e-04 7.243900568877903e-02   -0.171572573090887     0.045364639214749     3.908577322669729     0.045364639224539     1.977012221173589     0.242699761577846    -0.000000000000000     0.000000000000002     0.000000000000002    -0.017111709048689   -11.443726188143213     0.000320923084906
	// 0.030000  2.024566513093719e-01 2.851515745296326e-28 -4.185602694850953e-23 9.858301422097954e-23 2.628289935843225e-02 4.128716262753067e+00 7.710665349450496e-03 1.162212672883882e-01   -0.306947423752770     0.199145554757876     3.562871601770279     0.199147089634617     1.887557045964513     0.379477318142794     0.000000000000000    -0.000000000000000    -0.000000000000000    -0.007803549644695    -8.290129193920482     0.008475554107027

    initV_B_n999.set_barchiAB(2.019508440767480e-01);
    initV_B_n999.set_simchiAB(ZERO);
    initV_B_n999.set_DeltaS(ZERO);
    initV_B_n999.set_DeltaT(ZERO);
    initV_B_n999.set_simDeltaB(5.253586960997034e-02);
    initV_B_n999.set_mu(5.721566801475077e+00);
    initV_B_n999.set_mAF(2.439528519551680e-04);
    initV_B_n999.set_d(7.243900568877903e-02);

    initV_B_n97.set_barchiAB(2.024566513093719e-01);
    initV_B_n97.set_simchiAB(ZERO);
    initV_B_n97.set_DeltaS(ZERO);
    initV_B_n97.set_DeltaT(ZERO);
    initV_B_n97.set_simDeltaB(2.628289935843225e-02);
    initV_B_n97.set_mu(4.128716262753067e+00);
    initV_B_n97.set_mAF(7.710665349450496e-03);
    initV_B_n97.set_d(1.162212672883882e-01);

	// C
	// 0.001000  1.704967692013200e-01 2.413359112013626e-30 1.682188545380227e-01 -4.635888312848086e-13 9.210213641760333e-29 3.435732916258100e+00 1.995019207676190e-12 1.514704400911580e-02   -0.232912590267705     0.005034800590002     3.988345562770387     0.005034800590002     1.997084265315409     0.175034808817695     0.000000000000000     0.167728730517760     0.000000000000154    -0.000000000000000    -6.872981605930469     0.000000000001644
	// 0.030000  1.887184112277540e-01 4.022180390981821e-30 1.381803827054838e-01 -2.362239931623809e-12 -6.534899751061710e-29 3.557011413537177e+00 -9.298512036204839e-15 9.034273248871363e-02   -0.336694285705645     0.155819534716066     3.644551802712305     0.155819534716067     1.909070926579813     0.327768041177345     0.000000000000000     0.125901390722187     0.000000000000717     0.000000000000000    -7.176785467533806    -0.000000000000007

    initV_C_n999.set_barchiAB(1.704967692013200e-01);
    initV_C_n999.set_simchiAB(ZERO);
    initV_C_n999.set_DeltaS(1.682188545380227e-01);
    initV_C_n999.set_DeltaT(ZERO);
    initV_C_n999.set_simDeltaB(ZERO);
    initV_C_n999.set_mu(3.435732916258100e+00);
    initV_C_n999.set_mAF(ZERO);
    initV_C_n999.set_d(1.514704400911580e-02);

    initV_C_n97.set_barchiAB(1.887184112277540e-01);
    initV_C_n97.set_simchiAB(ZERO);
    initV_C_n97.set_DeltaS(1.381803827054838e-01);
    initV_C_n97.set_DeltaT(ZERO);
    initV_C_n97.set_simDeltaB(ZERO);
    initV_C_n97.set_mu(3.557011413537177e+00);
    initV_C_n97.set_mAF(ZERO);
    initV_C_n97.set_d(9.034273248871363e-02);

	// D
	// 0.001000  1.675985933924276e-01 2.053172127425808e-30 4.524268332861422e-02 1.130700608008279e-03 4.808200360334780e-27 5.683004924143524e+00 6.556539893009282e-01 8.362154625640089e-02   -0.206249363203051     0.072311449324424     1.898216238343181     0.072373755925106     1.377757684915305     0.151846042199670     0.000000000000000     0.021470099040148    -0.000178859521235    -0.000000000000000   -11.382637863648705     0.388872826373376
	// 0.030000  1.788807732537705e-01 3.852776015660254e-27 5.346195652101027e-02 2.347003866889418e-02 3.722191632085162e-27 3.361444393505849e+00 5.537974035700060e-01 1.138843887832251e-01   -0.323939198799920     0.209783114095174     2.090798137462619     0.217903743935764     1.445959244744685     0.303284010981888     0.000000000000000     0.027944539779809    -0.004089259427925    -0.000000000000000    -7.139273894109552     0.560016706500568

    initV_D_n999.set_barchiAB(1.675985933924276e-01);
    initV_D_n999.set_simchiAB(ZERO);
    initV_D_n999.set_DeltaS(4.524268332861422e-02);
    initV_D_n999.set_DeltaT(1.130700608008279e-03);
    initV_D_n999.set_simDeltaB(ZERO);
    initV_D_n999.set_mu(5.683004924143524e+00);
    initV_D_n999.set_mAF(6.556539893009282e-01);
    initV_D_n999.set_d(8.362154625640089e-02);

    initV_D_n97.set_barchiAB(1.788807732537705e-01);
    initV_D_n97.set_simchiAB(ZERO);
    initV_D_n97.set_DeltaS(5.346195652101027e-02);
    initV_D_n97.set_DeltaT(2.347003866889418e-02);
    initV_D_n97.set_simDeltaB(ZERO);
    initV_D_n97.set_mu(3.361444393505849e+00);
    initV_D_n97.set_mAF(5.537974035700060e-01);
    initV_D_n97.set_d(1.138843887832251e-01);

	// E
    //0.001000  1.693211838431733e-01 -9.010683092846660e-27 1.669113860227631e-01 3.921232154000801e-05 1.795614877775459e-27 3.371546654213554e+00 1.319009092615469e-01 1.496850257060662e-02   -0.233035552832709     0.004971366186468     3.853190565855477     0.004999079730841     1.962954550124754     0.168078063232467    -0.000000000000000     0.160785344489193    -0.000012591045619    -0.000000000000000    -6.744734991792987     0.100877447601264
    //0.030000  1.887184112277357e-01 7.153500092116778e-29 1.381803827054658e-01 -5.252053833144162e-14 1.059544162527236e-27 3.557011413536724e+00 -6.096904527941155e-13 9.034273248869337e-02   -0.336694285705647     0.155819534716035     3.644551802712361     0.155819534716035     1.909070926579827     0.327768041177300     0.000000000000000     0.125901390722172     0.000000000000016    -0.000000000000000    -7.176785467532885    -0.000000000000431

    initV_E_n999.set_barchiAB(1.693211838431733e-01);
    initV_E_n999.set_simchiAB(ZERO);
    initV_E_n999.set_DeltaS(1.669113860227631e-01);
    initV_E_n999.set_DeltaT(3.921232154000801e-05);
    initV_E_n999.set_simDeltaB(ZERO);
    initV_E_n999.set_mu(3.371546654213554e+00);
    initV_E_n999.set_mAF(1.319009092615469e-01);
    initV_E_n999.set_d(1.496850257060662e-02);

    initV_E_n97.set_barchiAB(1.887184112277357e-01);
    initV_E_n97.set_simchiAB(ZERO);
    initV_E_n97.set_DeltaS(1.381803827054658e-01);
    initV_E_n97.set_DeltaT(ZERO);
    initV_E_n97.set_simDeltaB(ZERO);
    initV_E_n97.set_mu(3.557011413536724e+00);
    initV_E_n97.set_mAF(ZERO);
    initV_E_n97.set_d(9.034273248869337e-02);

    // 5.000000 0.700000 1.907608990935007e-01 1.589558430705757e-28 2.310807918237033e-23 3.193343627839878e-27 4.028031263888447e-28 5.655897444477536e-01 5.018664661491306e-15 2.285326945938115e-01   -1.150943256154826     0.884440074793620     1.713196858107756     0.884440074793620     1.308891461545898
    initV_E_n3.set_barchiAB(1.907608990935007e-01);
    initV_E_n3.set_simchiAB(ZERO);
    initV_E_n3.set_DeltaS(ZERO);
    initV_E_n3.set_DeltaT(ZERO);
    initV_E_n3.set_simDeltaB(ZERO);
    initV_E_n3.set_mu(5.655897444477536e-01);
    initV_E_n3.set_mAF(ZERO);
    initV_E_n3.set_d(2.285326945938115e-01);

	Var initV = initV_E_n3;

	Res res(N_VAR);
	RootGSL root(par, function, res);

	initV.shake();
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      "); res.print_label();
	fflush(stdout);

	int LOOPS1 = 10;
	int LOOPS2 = 10;
	int LOOPS3 = 10;
	int LOOPS4 = 10;
	initV = old;
	for(int i=0; i<LOOPS1; i++) {
		par.set_L(128 + 16*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<=LOOPS2; i++) {
		par.set_L(288 + 32*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<LOOPS3; i++) {
		par.set_L(640 + 64*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<=LOOPS4; i++) {
		par.set_L(1280 + 128*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
}

void run_1D_J() {
    Par par(N_PAR);
    par.set_n(0.9);
    par.set_J(1./3.);
    par.set_t(-1.);
    par.set_U(100);
    par.set_beta(1500.);
    par.set_L(512);
	printf("#J-run 2D2t_SGA_tPrim0\n");
	printf("#"); par.print_label();
	printf("#"); par.print_verboten();

    double ZERO = 1e-7;

    // #sdelta  U           barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm
    //0.996000 30.000000 1.712339149661826e-01 -7.116046177010064e-25 1.667588557806268e-01 8.099181961153371e-05 2.475126730413155e-24 2.543659754840674e+00 6.868129192639835e-02 6.918257610033627e-03   -0.239724543090433     0.009881491279315     3.930261184553631     0.009941314187715     1.982488634154969
    Var initV_U30_J13_delta996(N_VAR);	// m+DS+DT
    initV_U30_J13_delta996.set_barchiAB(1.712339149661826e-01);
    initV_U30_J13_delta996.set_simchiAB(ZERO);
    initV_U30_J13_delta996.set_DeltaS(1.667588557806268e-01);
    initV_U30_J13_delta996.set_DeltaT(8.099181961153371e-05);
    initV_U30_J13_delta996.set_simDeltaB(ZERO);
    initV_U30_J13_delta996.set_mu(2.543659754840674e+00);
    initV_U30_J13_delta996.set_mAF(6.868129192639835e-02);
    initV_U30_J13_delta996.set_d(6.918257610033627e-03);

    // U=100
	//#sdelta          barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm                 lbarchiAB                 lsimchiAB             lDeltaS             lDeltaT              lsimDeltaB            ln                  lmAF
    //0.001000  1.686360066445106e-01 2.207533802489830e-28 1.676171218402892e-01 3.872342472544926e-05 1.204180861839948e-28 2.266511470421953e+00 1.319374026200395e-01 9.082992468738691e-04   -0.231943172431472     0.002083912304764     3.856563937898645     0.002149468211930     1.963813620967795     0.164672797768878     0.000000000000000     0.161606536865906    -0.000012444946946    -0.000000000000000    -4.534645832983689     0.100128507152207
    Var initV_U100_J13_delta999(N_VAR);	// m+DS+DT
    initV_U100_J13_delta999.set_barchiAB(1.686360066445106e-01);
    initV_U100_J13_delta999.set_simchiAB(ZERO);
    initV_U100_J13_delta999.set_DeltaS(1.676171218402892e-01);
    initV_U100_J13_delta999.set_DeltaT(3.872342472544926e-05);
    initV_U100_J13_delta999.set_simDeltaB(ZERO);
    initV_U100_J13_delta999.set_mu(2.266511470421953e+00);
    initV_U100_J13_delta999.set_mAF(1.319374026200395e-01);
    initV_U100_J13_delta999.set_d(9.082992468738691e-04);

    // U=100
	//#sdelta          barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm                 lbarchiAB                 lsimchiAB             lDeltaS             lDeltaT              lsimDeltaB            ln                  lmAF
    //0.010000  1.737023693201609e-01 5.933658230033263e-30 1.646250241932819e-01 -7.557268807420535e-15 1.781384152967369e-29 2.349572567029063e+00 2.584824515397958e-13 2.886053425307017e-03   -0.252892615073981     0.020978081913192     3.921052236165468     0.020978081913192     1.980164699252430     0.191252097825706     0.000000000000000     0.161375829810466     0.000000000000002    -0.000000000000000    -4.714968219169414     0.000000000000184
    Var initV_U100_J13_delta99(N_VAR);	// m+DS+DT
    initV_U100_J13_delta99.set_barchiAB(1.737023693201609e-01);
    initV_U100_J13_delta99.set_simchiAB(ZERO);
    initV_U100_J13_delta99.set_DeltaS(1.646250241932819e-01);
    initV_U100_J13_delta99.set_DeltaT(ZERO);
    initV_U100_J13_delta99.set_simDeltaB(ZERO);
    initV_U100_J13_delta99.set_mu(2.349572567029063e+00);
    initV_U100_J13_delta99.set_mAF(ZERO);
    initV_U100_J13_delta99.set_d(2.886053425307017e-03);

    // U=100
	//#sdelta          barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm                 lbarchiAB                 lsimchiAB             lDeltaS             lDeltaT              lsimDeltaB            ln                  lmAF
    //0.100000  1.906807046783698e-01 8.148002374013715e-27 1.272486700252829e-01 3.754395163827864e-13 -4.405400101762304e-25 2.165786305050962e+00 3.689082922803584e-14 9.157466575081923e-03   -0.459107195113225     0.192621982398796     3.304553147763666     0.192621982398796     1.817842993155258     0.350150613114471     0.000000000000000     0.105124998270197    -0.000000000000103     0.000000000000000    -4.540582591353240    -0.000000000000009
    Var initV_U100_J13_delta9(N_VAR);	// m+DS+DT
    initV_U100_J13_delta9.set_barchiAB(1.906807046783698e-01);
    initV_U100_J13_delta9.set_simchiAB(ZERO);
    initV_U100_J13_delta9.set_DeltaS(1.272486700252829e-01);
    initV_U100_J13_delta9.set_DeltaT(ZERO);
    initV_U100_J13_delta9.set_simDeltaB(ZERO);
    initV_U100_J13_delta9.set_mu(2.165786305050962e+00);
    initV_U100_J13_delta9.set_mAF(ZERO);
    initV_U100_J13_delta9.set_d(9.157466575081923e-03);

	// I'm choosing start point
	Var initV = initV_U100_J13_delta9;

	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#J           "); res.print_label();
	fflush(stdout);


	LOOPS = 100;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_J(1./3. - ((double)i)*1./3./LOOPS);
		printf("%f ", par.J());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 200;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_J(1./3. + ((double)i)*2./3./LOOPS);
		printf("%f ", par.J());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	LOOPS = 100;
	for(int i=0; i<=LOOPS; i++) {
		par.set_J(1. + ((double)i)*1./LOOPS);
		printf("%f ", par.J());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
};

void run_optimal_doping() {
    Par par(N_PAR);
    par.set_n(0.996);
    par.set_J(1./3.);
    par.set_t(-1.);
    par.set_U(100);
    par.set_beta(1500.);
    par.set_L(512);
	printf("#J-run 2D2t_SGA_tPrim0\n");
	printf("#"); par.print_label();
	printf("#"); par.print_verboten();

    double ZERO = 1e-7;

    // #sdelta  U           barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm
    //0.996000 30.000000 1.712339149661826e-01 -7.116046177010064e-25 1.667588557806268e-01 8.099181961153371e-05 2.475126730413155e-24 2.543659754840674e+00 6.868129192639835e-02 6.918257610033627e-03   -0.239724543090433     0.009881491279315     3.930261184553631     0.009941314187715     1.982488634154969
    Var initV_U30_J13_delta996(N_VAR);	// m+DS+DT
    initV_U30_J13_delta996.set_barchiAB(1.712339149661826e-01);
    initV_U30_J13_delta996.set_simchiAB(ZERO);
    initV_U30_J13_delta996.set_DeltaS(1.667588557806268e-01);
    initV_U30_J13_delta996.set_DeltaT(8.099181961153371e-05);
    initV_U30_J13_delta996.set_simDeltaB(ZERO);
    initV_U30_J13_delta996.set_mu(2.543659754840674e+00);
    initV_U30_J13_delta996.set_mAF(6.868129192639835e-02);
    initV_U30_J13_delta996.set_d(6.918257610033627e-03);

    //#J                            barchiAB             simchiAB                 DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm                 lbarchiAB                 lsimchiAB             lDeltaS             lDeltaT              lsimDeltaB            ln                  lmAF
    //0.333333 0.825990 28 1.940625553700410e-01 -4.582370558159427e-295 9.556682369377588e-02 2.735780860679661e-16 2.519572312785578e-266 1.867509767223852e+00 -1.776749157664378e-13 1.143644774390575e-02   -0.608634984514793     0.313043445707252     2.900290188945882     0.313043445707252     1.703023836869549     0.453752877052626    -0.000000000000000     0.069292880286945    -0.000000000000000    -0.000000000000000    -4.141900799395148     0.000000000000173
    Var initV_n999_max(N_VAR);	// m+DS+DT
    initV_n999_max.set_barchiAB(1.940608256504650e-01);
    initV_n999_max.set_simchiAB(ZERO);
    initV_n999_max.set_DeltaS(9.559213545597220e-02);
    initV_n999_max.set_DeltaT(ZERO);
    initV_n999_max.set_simDeltaB(ZERO);
    initV_n999_max.set_mu(1.867535760031845e+00);
    initV_n999_max.set_mAF(ZERO);
    initV_n999_max.set_d(1.143633442901475e-02);


	// I'm choosing start point
	Var initV = initV_U30_J13_delta996;
    //Var initV = initV_n999_max;

	Res res(N_VAR);
	RootGSL root(par, function, res);

	if(root.Solve(initV) != GSL_SUCCESS) { printf("#First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#J           "); res.print_label();
	fflush(stdout);

	// declarations
	//double max_DeltaS = res.get_V()->DeltaS()*glob.gBDelta;
	double current_n = par.n();
	int LOOPS = 200;

	for(int i=0; i<=LOOPS; i++) {
		par.set_J(1./3. + ((double)i)*2./3./LOOPS);
		//initV.shake();
		old = initV;

		double MAX_STEP;
		if(i==0) { MAX_STEP = 0.04; } else { MAX_STEP = 0.001; }
		double n_step = MAX_STEP;
		double MIN_STEP = 0.00001;

		double up_DeltaS, down_DeltaS;
		double max_DeltaS = 0;
		int j=0;

		// Main loop
		while(1) {
			j++;
			Var varup(N_VAR);
			Var vardown(N_VAR);
			initV.shake();

			double n_up = current_n + n_step;
			if(n_up <= 1.) {
				par.set_n(n_up);
				//fprintf(stdout, "#%.12f -> %.12f ( %.12f)\n", current_n, n_up, n_step); fflush(stdout);
				if(root.Solve(initV) != GSL_SUCCESS) { printf("#GSL FAILED\n"); initV = varup; fflush(stdout); break; }
				up_DeltaS = res.get_V()->DeltaS()*glob.gBDelta;
				varup.set(res.get_V());
			}
			else { up_DeltaS = 0;}

			double n_down = current_n - n_step;
			if(n_down >= 0.) {
				par.set_n(n_down);
				//fprintf(stdout, "#%.12f -> %.12f (-%.12f)\n", current_n, n_down, n_step); fflush(stdout);
				if(root.Solve(initV) != GSL_SUCCESS) { printf("#GSL FAILED\n"); initV = vardown; fflush(stdout); break; }
				down_DeltaS = res.get_V()->DeltaS()*glob.gBDelta;
				vardown.set(res.get_V());
			}
			else { down_DeltaS = 0;}

			//fprintf(stdout, "#%.12f <- %.12f -> %.12f | (-%.12f) | %.12f %.12f %.12f \n", n_down, current_n, n_up, n_step, down_DeltaS, max_DeltaS, up_DeltaS); fflush(stdout);

			if(up_DeltaS > max_DeltaS && up_DeltaS > down_DeltaS) {
				max_DeltaS = up_DeltaS;
				initV.set(res.get_V());
				current_n = n_up;
				initV = varup;
				continue;
			}
			else if (down_DeltaS > max_DeltaS && up_DeltaS < down_DeltaS) {
				max_DeltaS = down_DeltaS;
				initV.set(res.get_V());
				current_n = n_down;
				initV = vardown;
				continue;
			}

			if(n_step/2. <= MIN_STEP) {
				old = initV;
				par.set_n(current_n);
				initV.shake();
				if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
				initV.set(res.get_V());

				printf("%f %f %d", par.J(), par.n(), j);
				root.print();
				fflush(stdout);
				break;
			}
			n_step /= 2.;
		}
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);

/*
	for(int i=0; i<=LOOPS; i++) {
		par.set_J(1./3. - ((double)i)*1./3./LOOPS);
		initV.shake();
		old = initV;

		double temp_DeltaS;
		if(i==0) {
			MAX_STEP = 0.04;
			MIN_STEP = 0.00001;
		}
		else {
			MAX_STEP = 0.0001;
			MIN_STEP = 0.0000001;
		}
		double n_step = MAX_STEP;
		int j=0;

		// I'm choosing direction
		double PROOBE_RANGE = 0.00001;
		par.set_n(current_n + PROOBE_RANGE);
		if(root.Solve(initV) != GSL_SUCCESS) { printf("#GSL FAILED at beginning.\n"); return; }
		double temp_DeltaS_up = res.get_V()->DeltaS()*glob.gBDelta;
		initV = old;

		par.set_n(current_n - PROOBE_RANGE);
		if(root.Solve(initV) != GSL_SUCCESS) { printf("#GSL FAILED at beginning.\n"); return; }
		double temp_DeltaS_down = res.get_V()->DeltaS()*glob.gBDelta;
		initV = old;

		double mul; // case temp_DeltaS_up > temp_DeltaS_down // go with n up
		if (temp_DeltaS_up < temp_DeltaS_down) { mul = -1.; } // go with n down
		else if (temp_DeltaS_up > temp_DeltaS_down) { mul = 1.; }
		else {
			double r = temp_DeltaS_up - max_DeltaS;
			if(r < RES && r > -RES) {
				printf("%f %f %d", par.J(), par.n(), j);
				root.print();
				continue;
			}
			else { printf("#MIN_STEP to small or other exotic error.\n"); return; }
		}
		//printf("%%f %f %f %f\n", temp_DeltaS_down, max_DeltaS, temp_DeltaS_up, mul);

		// Main loop
		while(1) {
			j++;
			double n_new = current_n + n_step*mul;
			if(n_new >= 0. && n_new <= 1.) {
				par.set_n(n_new);
				fprintf(stdout, "#%f -> %f (%f)\n", current_n, n_new, n_step*mul); fflush(stdout);
				if(root.Solve(initV) != GSL_SUCCESS) { printf("#GSL FAILED\n"); initV = old; fflush(stdout); break; }
				temp_DeltaS = res.get_V()->DeltaS()*glob.gBDelta;
			}
			else { printf("#Out of range\n"); initV = old; fflush(stdout); break; }

			if(temp_DeltaS > max_DeltaS) {
				max_DeltaS = temp_DeltaS;
				initV.set(res.get_V());
				current_n = n_new;
				continue;
			}

			if(n_step <= MIN_STEP) {
				printf("%f %f %d", par.J(), par.n(), j);
				root.print();
				break;
			}
			n_step /= 2.;
		}
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
*/

	/*
	int j=0;
	while(1){
		j++;
		par.set_n(current_n - n_step);
		fprintf(stderr, "#%f -> %f (-%f)\n", current_n, current_n - n_step, n_step); fflush(stderr);
		if(root.Solve(initV) != GSL_SUCCESS) { printf("GSL FAILED\n"); initV = old; fflush(stdout); continue; }

		double temp_DeltaS = res.get_V()->DeltaS()*glob.gBDelta;
		if(temp_DeltaS > max_DeltaS) {
			max_DeltaS = temp_DeltaS;
			initV.set(res.get_V());
			current_n -= n_step;
			continue;
		}

		if(n_step < MIN_STEP) {
			printf("%f %f %d", par.J(), par.n(), j);
			root.print();
			break;
		}
		n_step /= 2.;
	}


	// another points
	int LOOPS;
	LOOPS = 100;
	double max_DeltaS = res.get_V()->DeltaS()*glob.gBDelta;
	double current_n = par.n();
	double MAX_STEP = 0.01;
	double MIN_STEP = 0.00001;
	double n_step = MAX_STEP;

	for(int i=0; i<=LOOPS; i++) {
		par.set_J(1./3. - ((double)i)*1./3./LOOPS);
		initV.shake();
		old = initV;

		int j=0;
		while(1) {
			j++;
			double temp_DeltaS_up, temp_DeltaS_down;

			//initV = old;
			double n_up = current_n + n_step;
			if(n_up <= 1.) {
				par.set_n(n_up);
				fprintf(stderr, "#%f -> %f (+%f)\n", current_n, n_up, n_step); fflush(stderr);
				if(root.Solve(initV) != GSL_SUCCESS) { printf("GSL FAILED\n"); initV = old; fflush(stdout); continue; }
				temp_DeltaS_up = res.get_V()->DeltaS()*glob.gBDelta;
			}
			else { printf("#Out of range\n"); initV = old; fflush(stdout); break; }

			if(temp_DeltaS_up > max_DeltaS) {
				max_DeltaS = temp_DeltaS_up;
				initV.set(res.get_V());
				current_n = n_up;
				continue;
			}

			initV = old;
			double n_down = current_n - n_step;
			if(n_down >= 0.) {
				par.set_n(n_down);
				fprintf(stderr, "#%f -> %f (-%f)\n", current_n, n_down, n_step); fflush(stderr);
				if(root.Solve(initV) != GSL_SUCCESS) { printf("GSL FAILED\n"); initV = old; fflush(stdout); continue; }
				temp_DeltaS_down = res.get_V()->DeltaS()*glob.gBDelta;
			}
			else { printf("#Out of range\n"); initV = old; fflush(stdout); break; }

			if(temp_DeltaS_down > max_DeltaS) {
				max_DeltaS = temp_DeltaS_down;
				initV.set(res.get_V());
				current_n = n_down;
				continue;
			}

			if(n_step <= MIN_STEP) {
				printf("%f %f %d", par.J(), par.n(), j);
				root.print();
				break;
			}
			n_step /= 2.;
		}
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
    */

};

} // end of namespace case2D2t_SGA_tPrim0
